home *** CD-ROM | disk | FTP | other *** search
/ Complete Linux / Complete Linux.iso / xwindows / demos / xfract_1.z / xfract_1 / xfractint-1.06 / fractals.c < prev    next >
C/C++ Source or Header  |  1992-09-28  |  88KB  |  3,380 lines
  1. /*
  2. FRACTALS.C, FRACTALP.C and CALCFRAC.C actually calculate the fractal
  3. images (well, SOMEBODY had to do it!).  The modules are set up so that
  4. all logic that is independent of any fractal-specific code is in
  5. CALCFRAC.C, the code that IS fractal-specific is in FRACTALS.C, and the
  6. struscture that ties (we hope!) everything together is in FRACTALP.C.
  7. Original author Tim Wegner, but just about ALL the authors have
  8. contributed SOME code to this routine at one time or another, or
  9. contributed to one of the many massive restructurings.
  10.  
  11. The Fractal-specific routines are divided into three categories:
  12.  
  13. 1. Routines that are called once-per-orbit to calculate the orbit
  14.    value. These have names like "XxxxFractal", and their function
  15.    pointers are stored in fractalspecific[fractype].orbitcalc. EVERY
  16.    new fractal type needs one of these. Return 0 to continue iterations,
  17.    1 if we're done. Results for integer fractals are left in 'lnew.x' and
  18.    'lnew.y', for floating point fractals in 'new.x' and 'new.y'.
  19.  
  20. 2. Routines that are called once per pixel to set various variables
  21.    prior to the orbit calculation. These have names like xxx_per_pixel
  22.    and are fairly generic - chances are one is right for your new type.
  23.    They are stored in fractalspecific[fractype].per_pixel.
  24.  
  25. 3. Routines that are called once per screen to set various variables.
  26.    These have names like XxxxSetup, and are stored in
  27.    fractalspecific[fractype].per_image.
  28.  
  29. 4. The main fractal routine. Usually this will be StandardFractal(),
  30.    but if you have written a stand-alone fractal routine independent
  31.    of the StandardFractal mechanisms, your routine name goes here,
  32.    stored in fractalspecific[fractype].calctype.per_image.
  33.  
  34. Adding a new fractal type should be simply a matter of adding an item
  35. to the 'fractalspecific' structure, writing (or re-using one of the existing)
  36. an appropriate setup, per_image, per_pixel, and orbit routines.
  37.  
  38. --------------------------------------------------------------------   */
  39.  
  40.  
  41. #include <stdio.h>
  42. #include <stdlib.h>
  43. #include <float.h>
  44. #include <limits.h>
  45. #include <string.h>
  46. #ifdef __TURBOC__
  47. #include <alloc.h>
  48. #else
  49. #include <malloc.h>
  50. #endif
  51. #include "fractint.h"
  52. #include "mpmath.h"
  53. #include "helpdefs.h"
  54. #include "fractype.h"
  55. #include "prototyp.h"
  56.  
  57. #define NEWTONDEGREELIMIT  100
  58.  
  59. extern struct complex  initorbit;
  60. extern struct lcomplex linitorbit;
  61. extern char useinitorbit;
  62. extern double fgLimit;
  63. extern int distest;
  64.  
  65. #define CMPLXsqr_old(out)    \
  66.    (out).y = (old.x+old.x) * old.y;\
  67.    (out).x = tempsqrx - tempsqry
  68.  
  69. #define CMPLXpwr(arg1,arg2,out)   (out)= ComplexPower((arg1), (arg2))
  70. #define CMPLXmult1(arg1,arg2,out)    Arg2->d = (arg1); Arg1->d = (arg2);\
  71.      dStkMul(); Arg1++; Arg2++; (out) = Arg2->d
  72. #define CMPLXmult(arg1,arg2,out)  \
  73.     {\
  74.        CMPLX TmP;\
  75.        TmP.x = (arg1).x*(arg2).x - (arg1).y*(arg2).y;\
  76.        TmP.y = (arg1).x*(arg2).y + (arg1).y*(arg2).x;\
  77.        (out) = TmP;\
  78.      }
  79. #define CMPLXadd(arg1,arg2,out)    \
  80.     (out).x = (arg1).x + (arg2).x; (out).y = (arg1).y + (arg2).y
  81. #define CMPLXsub(arg1,arg2,out)    \
  82.     (out).x = (arg1).x - (arg2).x; (out).y = (arg1).y - (arg2).y
  83. #define CMPLXtimesreal(arg,real,out)   \
  84.     (out).x = (arg).x*(real);\
  85.     (out).y = (arg).y*(real)
  86.  
  87. #define CMPLXrecip(arg,out)    \
  88.    { double denom; denom = sqr((arg).x) + sqr((arg).y);\
  89.      if(denom==0.0) {(out).x = 1.0e10;(out).y = 1.0e10;}else\
  90.     { (out).x =  (arg).x/denom;\
  91.      (out).y = -(arg).y/denom;}}
  92.  
  93. extern int xshift, yshift;
  94. extern int biomorph;
  95. extern int forcesymmetry;
  96. extern int symmetry;
  97. LCMPLX lcoefficient,lold,lnew,lparm, linit,ltmp,ltmp2,lparm2;
  98. long ltempsqrx,ltempsqry;
  99. extern int decomp[];
  100. extern double param[];
  101. extern int potflag;                   /* potential enabled? */
  102. extern double f_radius,f_xcenter,f_ycenter;    /* inversion radius, center */
  103. extern double xxmin,xxmax,yymin,yymax;           /* corners */
  104. extern int overflow;
  105. extern int integerfractal;    /* TRUE if fractal uses integer math */
  106.  
  107. int maxcolor;
  108. int root, degree,basin;
  109. double floatmin,floatmax;
  110. double roverd, d1overd, threshold;
  111. CMPLX tmp2;
  112. extern CMPLX init,tmp,old,new,saved;
  113. CMPLX    staticroots[16]; /* roots array for degree 16 or less */
  114. CMPLX    *roots = staticroots;
  115. struct MPC    *MPCroots;
  116. extern int color, row, col;
  117. extern int invert;
  118. extern double far *dx0, far *dy0;
  119. extern double far *dx1, far *dy1;
  120. long FgHalf;
  121.  
  122. CMPLX one;
  123. CMPLX pwr;
  124. CMPLX Coefficient;
  125.  
  126. extern int    colors;             /* maximum colors available */
  127. extern int    inside;             /* "inside" color to use    */
  128. extern int    outside;            /* "outside" color to use   */
  129. extern int    finattract;
  130. extern int    fractype;            /* fractal type */
  131. extern int    debugflag;            /* for debugging purposes */
  132.  
  133. extern double    param[];        /* parameters */
  134. extern long    far *lx0, far *ly0;        /* X, Y points */
  135. extern long    far *lx1, far *ly1;        /* X, Y points */
  136. extern long    delx,dely;            /* X, Y increments */
  137. extern long    delmin;             /* min(max(dx),max(dy) */
  138. extern double    ddelmin;            /* min(max(dx),max(dy) */
  139. extern long    fudge;                /* fudge factor (2**n) */
  140. extern int    bitshift;            /* bit shift for fudge */
  141. int    bitshiftless1;            /* bit shift less 1 */
  142.  
  143. #ifndef sqr
  144. #define sqr(x) ((x)*(x))
  145. #endif
  146.  
  147. #ifndef lsqr
  148. #define lsqr(x) (multiply((x),(x),bitshift))
  149. #endif
  150.  
  151. #define modulus(z)     (sqr((z).x)+sqr((z).y))
  152. #define conjugate(pz)    ((pz)->y = 0.0 - (pz)->y)
  153. #define distance(z1,z2)  (sqr((z1).x-(z2).x)+sqr((z1).y-(z2).y))
  154. #define pMPsqr(z) (*pMPmul((z),(z)))
  155. #define MPdistance(z1,z2)  (*pMPadd(pMPsqr(*pMPsub((z1).x,(z2).x)),pMPsqr(*pMPsub((z1).y,(z2).y))))
  156.  
  157. double twopi = PI*2.0;
  158. static int c_exp;
  159.  
  160.  
  161. /* These are local but I don't want to pass them as parameters */
  162. CMPLX lambda;
  163. extern double magnitude, rqlim, rqlim2;
  164. CMPLX parm,parm2;
  165. CMPLX *floatparm;
  166. LCMPLX *longparm; /* used here and in jb.c */
  167. extern int (*calctype)();
  168. extern unsigned long lm;        /* magnitude limit (CALCMAND) */
  169.  
  170. /* -------------------------------------------------------------------- */
  171. /*        These variables are external for speed's sake only      */
  172. /* -------------------------------------------------------------------- */
  173.  
  174. double sinx,cosx,sinhx,coshx;
  175. double siny,cosy,sinhy,coshy;
  176. double tmpexp;
  177. double tempsqrx,tempsqry;
  178.  
  179. double foldxinitx,foldyinity,foldxinity,foldyinitx;
  180. long oldxinitx,oldyinity,oldxinity,oldyinitx;
  181. long longtmp;
  182. extern long lmagnitud, llimit, llimit2, l16triglim;
  183. extern periodicitycheck;
  184. extern char floatflag;
  185.  
  186. /* These are for quaternions */
  187. double qc,qci,qcj,qck;
  188.  
  189. /* temporary variables for trig use */
  190. long lcosx, lcoshx, lsinx, lsinhx;
  191. long lcosy, lcoshy, lsiny, lsinhy;
  192.  
  193. /*
  194. **  details of finite attractors (required for Magnet Fractals)
  195. **  (can also be used in "coloring in" the lakes of Julia types)
  196. */
  197. extern          int      attractors; /* number of finite attractors   */
  198. extern CMPLX  attr[];       /* finite attractor values (f.p) */
  199. extern LCMPLX lattr[];      /* finite attractor values (int) */
  200. extern int    attrperiod[]; /* finite attractor period */
  201.  
  202. /*
  203. **  pre-calculated values for fractal types Magnet2M & Magnet2J
  204. */
  205. CMPLX    T_Cm1;          /* 3 * (floatparm - 1)            */
  206. CMPLX    T_Cm2;          /* 3 * (floatparm - 2)            */
  207. CMPLX    T_Cm1Cm2;     /* (floatparm - 1) * (floatparm - 2) */
  208.  
  209. void FloatPreCalcMagnet2() /* precalculation for Magnet2 (M & J) for speed */
  210.   {
  211.     T_Cm1.x = floatparm->x - 1.0;   T_Cm1.y = floatparm->y;
  212.     T_Cm2.x = floatparm->x - 2.0;   T_Cm2.y = floatparm->y;
  213.     T_Cm1Cm2.x = (T_Cm1.x * T_Cm2.x) - (T_Cm1.y * T_Cm2.y);
  214.     T_Cm1Cm2.y = (T_Cm1.x * T_Cm2.y) + (T_Cm1.y * T_Cm2.x);
  215.     T_Cm1.x += T_Cm1.x + T_Cm1.x;   T_Cm1.y += T_Cm1.y + T_Cm1.y;
  216.     T_Cm2.x += T_Cm2.x + T_Cm2.x;   T_Cm2.y += T_Cm2.y + T_Cm2.y;
  217.   }
  218.  
  219.  
  220. /* -------------------------------------------------------------------- */
  221. /*        Bailout Routines Macros                                                 */
  222. /* -------------------------------------------------------------------- */
  223.  
  224. static int near floatbailout()
  225. {
  226.    if ( ( magnitude = ( tempsqrx=sqr(new.x) )
  227.             + ( tempsqry=sqr(new.y) ) ) >= rqlim ) return(1);
  228.    old = new;
  229.    return(0);
  230. }
  231.  
  232. /* longbailout() is equivalent to next */
  233. #define LONGBAILOUT()    \
  234.    ltempsqrx = lsqr(lnew.x); ltempsqry = lsqr(lnew.y);\
  235.    lmagnitud = ltempsqrx + ltempsqry;\
  236.    if (lmagnitud >= llimit || lmagnitud < 0 || labs(lnew.x) > llimit2\
  237.      || labs(lnew.y) > llimit2 || overflow) \
  238.            { overflow=0;return(1);}\
  239.    lold = lnew;
  240.  
  241. #define FLOATTRIGBAILOUT()  \
  242.    if (fabs(old.y) >= rqlim2) return(1);
  243.  
  244. #define LONGTRIGBAILOUT()  \
  245.    if(labs(lold.y) >= llimit2 || overflow) { overflow=0;return(1);}
  246.  
  247. #define LONGXYTRIGBAILOUT()  \
  248.    if(labs(lold.x) >= llimit2 || labs(lold.y) >= llimit2 || overflow)\
  249.     { overflow=0;return(1);}
  250.  
  251. #define FLOATXYTRIGBAILOUT()  \
  252.    if (fabs(old.x) >= rqlim2 || fabs(old.y) >= rqlim2) return(1);
  253.  
  254. #define FLOATHTRIGBAILOUT()  \
  255.    if (fabs(old.x) >= rqlim2) return(1);
  256.  
  257. #define LONGHTRIGBAILOUT()  \
  258.    if(labs(lold.x) >= llimit2 || overflow) { overflow=0;return(1);}
  259.  
  260. #define TRIG16CHECK(X)    \
  261.       if(labs((X)) > l16triglim || overflow) { overflow=0;return(1);}
  262.  
  263. #define FLOATEXPBAILOUT()  \
  264.    if (fabs(old.y) >= 1.0e8) return(1);\
  265.    if (fabs(old.x) >= 6.4e2) return(1);
  266.  
  267. #define LONGEXPBAILOUT()  \
  268.    if (labs(lold.y) >= (1000L<<bitshift)) return(1);\
  269.    if (labs(lold.x) >=      (8L<<bitshift)) return(1);
  270.  
  271. #if 0
  272. /* this define uses usual trig instead of fast trig */
  273. #define FPUsincos(px,psinx,pcosx) \
  274.    *(psinx) = sin(*(px));\
  275.    *(pcosx) = cos(*(px));
  276.  
  277. #define FPUsinhcosh(px,psinhx,pcoshx) \
  278.    *(psinhx) = sinh(*(px));\
  279.    *(pcoshx) = cosh(*(px));
  280. #endif
  281.  
  282. #define LTRIGARG(X)    \
  283.    if(labs((X)) > l16triglim)\
  284.    {\
  285.       double tmp;\
  286.       tmp = (X);\
  287.       tmp /= fudge;\
  288.       tmp = fmod(tmp,twopi);\
  289.       tmp *= fudge;\
  290.       (X) = tmp;\
  291.    }\
  292.  
  293. static int near Halleybailout()
  294. {
  295.    if ( fabs(modulus(new)-modulus(old)) < parm2.x)
  296.       return(1);
  297.    old = new;
  298.    return(0);
  299. }
  300.  
  301. #ifndef XFRACT
  302. #define MPCmod(m) (*pMPadd(*pMPmul((m).x, (m).x), *pMPmul((m).y, (m).y)))
  303. struct MPC mpcold, mpcnew, mpctmp, mpctmp1;
  304. struct MP mptmpparm2x;
  305.  
  306. static int near MPCHalleybailout()
  307. {
  308. static struct MP mptmpbailout;
  309.    mptmpbailout = *MPabs(*pMPsub(MPCmod(mpcnew), MPCmod(mpcold)));
  310.    if (pMPcmp(mptmpbailout, mptmpparm2x) < 0)
  311.       return(1);
  312.    mpcold = mpcnew;
  313.    return(0);
  314. }
  315. #endif
  316.  
  317. /* -------------------------------------------------------------------- */
  318. /*        Fractal (once per iteration) routines            */
  319. /* -------------------------------------------------------------------- */
  320. static double xt, yt, t2;
  321.  
  322. /* Raise complex number (base) to the (exp) power, storing the result
  323. ** in complex (result).
  324. */
  325. void cpower(CMPLX *base, int exp, CMPLX *result)
  326. {
  327.     if (exp<0) {
  328.     cpower(base,-exp,result);
  329.     CMPLXrecip(*result,*result);
  330.     return;
  331.     }
  332.  
  333.     xt = base->x;   yt = base->y;
  334.  
  335.     if (exp & 1)
  336.     {
  337.        result->x = xt;
  338.        result->y = yt;
  339.     }
  340.     else
  341.     {
  342.        result->x = 1.0;
  343.        result->y = 0.0;
  344.     }
  345.  
  346.     exp >>= 1;
  347.     while (exp)
  348.     {
  349.     t2 = xt * xt - yt * yt;
  350.     yt = 2 * xt * yt;
  351.     xt = t2;
  352.  
  353.     if (exp & 1)
  354.     {
  355.         t2 = xt * result->x - yt * result->y;
  356.         result->y = result->y * xt + yt * result->x;
  357.         result->x = t2;
  358.     }
  359.     exp >>= 1;
  360.     }
  361. }
  362. /* long version */
  363. static long lxt, lyt, lt2;
  364.  
  365. lcpower(LCMPLX *base, int exp, LCMPLX *result, int bitshift)
  366. {
  367.     static long maxarg;
  368.     maxarg = 64L<<bitshift;
  369.  
  370.     overflow = 0;
  371.     lxt = base->x;   lyt = base->y;
  372.  
  373.     if (exp & 1)
  374.     {
  375.        result->x = lxt;
  376.        result->y = lyt;
  377.     }
  378.     else
  379.     {
  380.        result->x = 1L<<bitshift;
  381.        result->y = 0L;
  382.     }
  383.  
  384.     exp >>= 1;
  385.     while (exp)
  386.     {
  387.     /*
  388.     if(labs(lxt) >= maxarg || labs(lyt) >= maxarg)
  389.        return(-1);
  390.     */
  391.     lt2 = multiply(lxt, lxt, bitshift) - multiply(lyt,lyt,bitshift);
  392.     lyt = multiply(lxt,lyt,bitshiftless1);
  393.     if(overflow)
  394.        return(overflow);
  395.     lxt = lt2;
  396.  
  397.     if (exp & 1)
  398.     {
  399.         lt2 = multiply(lxt,result->x, bitshift) - multiply(lyt,result->y,bitshift);
  400.         result->y = multiply(result->y,lxt,bitshift) + multiply(lyt,result->x,bitshift);
  401.         result->x = lt2;
  402.     }
  403.     exp >>= 1;
  404.     }
  405.     if(result->x == 0 && result->y == 0)
  406.        overflow = 1;
  407.     return(overflow);
  408. }
  409. #if 0
  410. z_to_the_z(CMPLX *z, CMPLX *out)
  411. {
  412.     static CMPLX tmp1,tmp2;
  413.     /* raises complex z to the z power */
  414.     int errno_xxx;
  415.     errno_xxx = 0;
  416.  
  417.     if(fabs(z->x) < DBL_EPSILON) return(-1);
  418.  
  419.     /* log(x + iy) = 1/2(log(x*x + y*y) + i(arc_tan(y/x)) */
  420.     tmp1.x = .5*log(sqr(z->x)+sqr(z->y));
  421.  
  422.     /* the fabs in next line added to prevent discontinuity in image */
  423.     tmp1.y = atan(fabs(z->y/z->x));
  424.  
  425.     /* log(z)*z */
  426.     tmp2.x = tmp1.x * z->x - tmp1.y * z->y;
  427.     tmp2.y = tmp1.x * z->y + tmp1.y * z->x;
  428.  
  429.     /* z*z = e**(log(z)*z) */
  430.     /* e**(x + iy) =  e**x * (cos(y) + isin(y)) */
  431.  
  432.     tmpexp = exp(tmp2.x);
  433.  
  434.     FPUsincos(&tmp2.y,&siny,&cosy);
  435.     out->x = tmpexp*cosy;
  436.     out->y = tmpexp*siny;
  437.     return(errno_xxx);
  438. }
  439. #endif
  440.  
  441. int complex_div(CMPLX arg1,CMPLX arg2,CMPLX *pz);
  442. int complex_mult(CMPLX arg1,CMPLX arg2,CMPLX *pz);
  443.  
  444. #ifdef XFRACT /* fractint uses the NewtonFractal2 code in newton.asm */
  445.  
  446. /* Distance of complex z from unit circle */
  447. #define DIST1(z) (((z).x-1.0)*((z).x-1.0)+((z).y)*((z).y))
  448. #define LDIST1(z) (lsqr((((z).x)-fudge)) + lsqr(((z).y)))
  449.  
  450.  
  451. int NewtonFractal2()
  452. {
  453.     static char start=1;
  454.     if(start)
  455.     {
  456.        start = 0;
  457.     }
  458.     cpower(&old, degree-1, &tmp);
  459.     complex_mult(tmp, old, &new);
  460.  
  461.     if (DIST1(new) < threshold)
  462.     {
  463.        if(fractype==NEWTBASIN || fractype==MPNEWTBASIN)
  464.        {
  465.       int tmpcolor;
  466.       int i;
  467.       tmpcolor = -1;
  468.       /* this code determines which degree-th root of root the
  469.          Newton formula converges to. The roots of a 1 are
  470.          distributed on a circle of radius 1 about the origin. */
  471.       for(i=0;i<degree;i++)
  472.          /* color in alternating shades with iteration according to
  473.         which root of 1 it converged to */
  474.           if(distance(roots[i],old) < threshold)
  475.           {
  476.           if (basin==2) {
  477.                   tmpcolor = 1+(i&7)+((color&1)<<3);
  478.           } else {
  479.               tmpcolor = 1+i;
  480.           }
  481.           break;
  482.           }
  483.        if(tmpcolor == -1)
  484.           color = maxcolor;
  485.        else
  486.           color = tmpcolor;
  487.        }
  488.        return(1);
  489.     }
  490.     new.x = d1overd * new.x + roverd;
  491.     new.y *= d1overd;
  492.  
  493.     /* Watch for divide underflow */
  494.     if ((t2 = tmp.x * tmp.x + tmp.y * tmp.y) < FLT_MIN)
  495.       return(1);
  496.     else
  497.     {
  498.     t2 = 1.0 / t2;
  499.     old.x = t2 * (new.x * tmp.x + new.y * tmp.y);
  500.     old.y = t2 * (new.y * tmp.x - new.x * tmp.y);
  501.     }
  502.     return(0);
  503. }
  504.  
  505.  
  506.  
  507. #endif /* newton code only used by xfractint */
  508.  
  509. complex_mult(arg1,arg2,pz)
  510. CMPLX arg1,arg2,*pz;
  511. {
  512.    pz->x = arg1.x*arg2.x - arg1.y*arg2.y;
  513.    pz->y = arg1.x*arg2.y+arg1.y*arg2.x;
  514.    return(0);
  515. }
  516.  
  517.  
  518. complex_div(numerator,denominator,pout)
  519. CMPLX numerator,denominator,*pout;
  520. {
  521.    double mod;
  522.    if((mod = modulus(denominator)) < FLT_MIN)
  523.       return(1);
  524.    conjugate(&denominator);
  525.    complex_mult(numerator,denominator,pout);
  526.    pout->x = pout->x/mod;
  527.    pout->y = pout->y/mod;
  528.    return(0);
  529. }
  530.  
  531. #ifndef XFRACT
  532. struct MP mproverd, mpd1overd, mpthreshold,sqrmpthreshold;
  533. struct MP mpt2;
  534. struct MP mpone;
  535. extern struct MPC MPCone;
  536. extern int MPOverflow;
  537. #endif
  538.  
  539. int MPCNewtonFractal()
  540. {
  541. #ifndef XFRACT
  542.     MPOverflow = 0;
  543.     mpctmp   = MPCpow(mpcold,degree-1);
  544.  
  545.     mpcnew.x = *pMPsub(*pMPmul(mpctmp.x,mpcold.x),*pMPmul(mpctmp.y,mpcold.y));
  546.     mpcnew.y = *pMPadd(*pMPmul(mpctmp.x,mpcold.y),*pMPmul(mpctmp.y,mpcold.x));
  547.     mpctmp1.x = *pMPsub(mpcnew.x, MPCone.x);
  548.     mpctmp1.y = *pMPsub(mpcnew.y, MPCone.y);
  549.     if(pMPcmp(MPCmod(mpctmp1),mpthreshold)< 0)
  550.     {
  551.       if(fractype==MPNEWTBASIN)
  552.       {
  553.      int tmpcolor;
  554.      int i;
  555.      tmpcolor = -1;
  556.      for(i=0;i<degree;i++)
  557.          if(pMPcmp(MPdistance(MPCroots[i],mpcold),mpthreshold) < 0)
  558.          {
  559.         if(basin==2)
  560.            tmpcolor = 1+(i&7) + ((color&1)<<3);
  561.         else
  562.            tmpcolor = 1+i;
  563.             break;
  564.          }
  565.       if(tmpcolor == -1)
  566.          color = maxcolor;
  567.       else
  568.          color = tmpcolor;
  569.       }
  570.        return(1);
  571.     }
  572.  
  573.     mpcnew.x = *pMPadd(*pMPmul(mpd1overd,mpcnew.x),mproverd);
  574.     mpcnew.y = *pMPmul(mpcnew.y,mpd1overd);
  575.     mpt2 = MPCmod(mpctmp);
  576.     mpt2 = *pMPdiv(mpone,mpt2);
  577.     mpcold.x = *pMPmul(mpt2,(*pMPadd(*pMPmul(mpcnew.x,mpctmp.x),*pMPmul(mpcnew.y,mpctmp.y))));
  578.     mpcold.y = *pMPmul(mpt2,(*pMPsub(*pMPmul(mpcnew.y,mpctmp.x),*pMPmul(mpcnew.x,mpctmp.y))));
  579.     new.x = *pMP2d(mpcold.x);
  580.     new.y = *pMP2d(mpcold.y);
  581.     return(MPOverflow);
  582. #endif
  583. }
  584.  
  585.  
  586. Barnsley1Fractal()
  587. {
  588. #ifndef XFRACT
  589.    /* Barnsley's Mandelbrot type M1 from "Fractals
  590.    Everywhere" by Michael Barnsley, p. 322 */
  591.  
  592.    /* calculate intermediate products */
  593.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  594.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  595.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  596.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  597.    /* orbit calculation */
  598.    if(lold.x >= 0)
  599.    {
  600.       lnew.x = (oldxinitx - longparm->x - oldyinity);
  601.       lnew.y = (oldyinitx - longparm->y + oldxinity);
  602.    }
  603.    else
  604.    {
  605.       lnew.x = (oldxinitx + longparm->x - oldyinity);
  606.       lnew.y = (oldyinitx + longparm->y + oldxinity);
  607.    }
  608.    return(longbailout());
  609. #endif
  610. }
  611.  
  612. Barnsley1FPFractal()
  613. {
  614.    /* Barnsley's Mandelbrot type M1 from "Fractals
  615.    Everywhere" by Michael Barnsley, p. 322 */
  616.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  617.  
  618.    /* calculate intermediate products */
  619.    foldxinitx = old.x * floatparm->x;
  620.    foldyinity = old.y * floatparm->y;
  621.    foldxinity = old.x * floatparm->y;
  622.    foldyinitx = old.y * floatparm->x;
  623.    /* orbit calculation */
  624.    if(old.x >= 0)
  625.    {
  626.       new.x = (foldxinitx - floatparm->x - foldyinity);
  627.       new.y = (foldyinitx - floatparm->y + foldxinity);
  628.    }
  629.    else
  630.    {
  631.       new.x = (foldxinitx + floatparm->x - foldyinity);
  632.       new.y = (foldyinitx + floatparm->y + foldxinity);
  633.    }
  634.    return(floatbailout());
  635. }
  636.  
  637. Barnsley2Fractal()
  638. {
  639. #ifndef XFRACT
  640.    /* An unnamed Mandelbrot/Julia function from "Fractals
  641.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  642.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  643.  
  644.    /* calculate intermediate products */
  645.    oldxinitx   = multiply(lold.x, longparm->x, bitshift);
  646.    oldyinity   = multiply(lold.y, longparm->y, bitshift);
  647.    oldxinity   = multiply(lold.x, longparm->y, bitshift);
  648.    oldyinitx   = multiply(lold.y, longparm->x, bitshift);
  649.  
  650.    /* orbit calculation */
  651.    if(oldxinity + oldyinitx >= 0)
  652.    {
  653.       lnew.x = oldxinitx - longparm->x - oldyinity;
  654.       lnew.y = oldyinitx - longparm->y + oldxinity;
  655.    }
  656.    else
  657.    {
  658.       lnew.x = oldxinitx + longparm->x - oldyinity;
  659.       lnew.y = oldyinitx + longparm->y + oldxinity;
  660.    }
  661.    return(longbailout());
  662. #endif
  663. }
  664.  
  665. Barnsley2FPFractal()
  666. {
  667.    /* An unnamed Mandelbrot/Julia function from "Fractals
  668.    Everywhere" by Michael Barnsley, p. 331, example 4.2 */
  669.  
  670.    /* calculate intermediate products */
  671.    foldxinitx = old.x * floatparm->x;
  672.    foldyinity = old.y * floatparm->y;
  673.    foldxinity = old.x * floatparm->y;
  674.    foldyinitx = old.y * floatparm->x;
  675.  
  676.    /* orbit calculation */
  677.    if(foldxinity + foldyinitx >= 0)
  678.    {
  679.       new.x = foldxinitx - floatparm->x - foldyinity;
  680.       new.y = foldyinitx - floatparm->y + foldxinity;
  681.    }
  682.    else
  683.    {
  684.       new.x = foldxinitx + floatparm->x - foldyinity;
  685.       new.y = foldyinitx + floatparm->y + foldxinity;
  686.    }
  687.    return(floatbailout());
  688. }
  689.  
  690. JuliaFractal()
  691. {
  692. #ifndef XFRACT
  693.    /* used for C prototype of fast integer math routines for classic
  694.       Mandelbrot and Julia */
  695.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  696.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  697.    return(longbailout());
  698. #else
  699.    fprintf(stderr,"JuliaFractal called\n");
  700.    exit(-1);
  701. #endif
  702. }
  703.  
  704. JuliafpFractal()
  705. {
  706.    /* floating point version of classical Mandelbrot/Julia */
  707.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  708.    new.x = tempsqrx - tempsqry + floatparm->x;
  709.    new.y = 2.0 * old.x * old.y + floatparm->y;
  710.    return(floatbailout());
  711. }
  712.  
  713. static double f(double x, double y)
  714. {
  715.    return(x - x*y);
  716. }
  717.  
  718. static double g(double x, double y)
  719. {
  720.    return(-y + x*y);
  721. }
  722.  
  723. LambdaFPFractal()
  724. {
  725.    /* variation of classical Mandelbrot/Julia */
  726.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  727.  
  728.    tempsqrx = old.x - tempsqrx + tempsqry;
  729.    tempsqry = -(old.y * old.x);
  730.    tempsqry += tempsqry + old.y;
  731.  
  732.    new.x = floatparm->x * tempsqrx - floatparm->y * tempsqry;
  733.    new.y = floatparm->x * tempsqry + floatparm->y * tempsqrx;
  734.    return(floatbailout());
  735. }
  736.  
  737. LambdaFractal()
  738. {
  739. #ifndef XFRACT
  740.    /* variation of classical Mandelbrot/Julia */
  741.  
  742.    /* in complex math) temp = Z * (1-Z) */
  743.    ltempsqrx = lold.x - ltempsqrx + ltempsqry;
  744.    ltempsqry = lold.y
  745.          - multiply(lold.y, lold.x, bitshiftless1);
  746.    /* (in complex math) Z = Lambda * Z */
  747.    lnew.x = multiply(longparm->x, ltempsqrx, bitshift)
  748.     - multiply(longparm->y, ltempsqry, bitshift);
  749.    lnew.y = multiply(longparm->x, ltempsqry, bitshift)
  750.     + multiply(longparm->y, ltempsqrx, bitshift);
  751.    return(longbailout());
  752. #endif
  753. }
  754.  
  755. SierpinskiFractal()
  756. {
  757. #ifndef XFRACT
  758.    /* following code translated from basic - see "Fractals
  759.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
  760.    lnew.x = (lold.x << 1);        /* new.x = 2 * old.x  */
  761.    lnew.y = (lold.y << 1);        /* new.y = 2 * old.y  */
  762.    if(lold.y > ltmp.y)    /* if old.y > .5 */
  763.       lnew.y = lnew.y - ltmp.x; /* new.y = 2 * old.y - 1 */
  764.    else if(lold.x > ltmp.y)    /* if old.x > .5 */
  765.       lnew.x = lnew.x - ltmp.x; /* new.x = 2 * old.x - 1 */
  766.    /* end barnsley code */
  767.    return(longbailout());
  768. #endif
  769. }
  770.  
  771. SierpinskiFPFractal()
  772. {
  773.    /* following code translated from basic - see "Fractals
  774.    Everywhere" by Michael Barnsley, p. 251, Program 7.1.1 */
  775.  
  776.    new.x = old.x + old.x;
  777.    new.y = old.y + old.y;
  778.    if(old.y > .5)
  779.       new.y = new.y - 1;
  780.    else if (old.x > .5)
  781.       new.x = new.x - 1;
  782.  
  783.    /* end barnsley code */
  784.    return(floatbailout());
  785. }
  786.  
  787. LambdaexponentFractal()
  788. {
  789.    /* found this in  "Science of Fractal Images" */
  790.    FLOATEXPBAILOUT();
  791.    FPUsincos  (&old.y,&siny,&cosy);
  792.  
  793.    if (old.x >= rqlim && cosy >= 0.0) return(1);
  794.    tmpexp = exp(old.x);
  795.    tmp.x = tmpexp*cosy;
  796.    tmp.y = tmpexp*siny;
  797.  
  798.    /*multiply by lamda */
  799.    new.x = floatparm->x*tmp.x - floatparm->y*tmp.y;
  800.    new.y = floatparm->y*tmp.x + floatparm->x*tmp.y;
  801.    old = new;
  802.    return(0);
  803. }
  804.  
  805. LongLambdaexponentFractal()
  806. {
  807. #ifndef XFRACT
  808.    /* found this in  "Science of Fractal Images" */
  809.    LONGEXPBAILOUT();
  810.  
  811.    SinCos086  (lold.y, &lsiny,    &lcosy);
  812.  
  813.    if (lold.x >= llimit && lcosy >= 0L) return(1);
  814.    longtmp = Exp086(lold.x);
  815.  
  816.    ltmp.x = multiply(longtmp,       lcosy,   bitshift);
  817.    ltmp.y = multiply(longtmp,       lsiny,   bitshift);
  818.  
  819.    lnew.x  = multiply(longparm->x, ltmp.x, bitshift)
  820.        - multiply(longparm->y, ltmp.y, bitshift);
  821.    lnew.y  = multiply(longparm->x, ltmp.y, bitshift)
  822.        + multiply(longparm->y, ltmp.x, bitshift);
  823.    lold = lnew;
  824.    return(0);
  825. #endif
  826. }
  827.  
  828. FloatTrigPlusExponentFractal()
  829. {
  830.    /* another Scientific American biomorph type */
  831.    /* z(n+1) = e**z(n) + trig(z(n)) + C */
  832.  
  833.    if (fabs(old.x) >= 6.4e2) return(1); /* DOMAIN errors */
  834.    tmpexp = exp(old.x);
  835.    FPUsincos  (&old.y,&siny,&cosy);
  836.    CMPLXtrig0(old,new);
  837.  
  838.    /*new =   trig(old) + e**old + C  */
  839.    new.x += tmpexp*cosy + floatparm->x;
  840.    new.y += tmpexp*siny + floatparm->y;
  841.    return(floatbailout());
  842. }
  843.  
  844.  
  845. LongTrigPlusExponentFractal()
  846. {
  847. #ifndef XFRACT
  848.    /* calculate exp(z) */
  849.  
  850.    /* domain check for fast transcendental functions */
  851.    TRIG16CHECK(lold.x);
  852.    TRIG16CHECK(lold.y);
  853.  
  854.    longtmp = Exp086(lold.x);
  855.    SinCos086  (lold.y, &lsiny,    &lcosy);
  856.    LCMPLXtrig0(lold,lnew);
  857.    lnew.x += multiply(longtmp,      lcosy,   bitshift) + longparm->x;
  858.    lnew.y += multiply(longtmp,      lsiny,   bitshift) + longparm->y;
  859.    return(longbailout());
  860. #endif
  861. }
  862.  
  863. MarksLambdaFractal()
  864. {
  865.    /* Mark Peterson's variation of "lambda" function */
  866.  
  867.    /* Z1 = (C^(exp-1) * Z**2) + C */
  868.    ltmp.x = ltempsqrx - ltempsqry;
  869.    ltmp.y = multiply(lold.x ,lold.y ,bitshiftless1);
  870.  
  871.    lnew.x = multiply(lcoefficient.x, ltmp.x, bitshift)
  872.     - multiply(lcoefficient.y, ltmp.y, bitshift) + longparm->x;
  873.    lnew.y = multiply(lcoefficient.x, ltmp.y, bitshift)
  874.     + multiply(lcoefficient.y, ltmp.x, bitshift) + longparm->y;
  875.  
  876.    return(longbailout());
  877. }
  878.  
  879.  
  880. long XXOne, FgOne, FgTwo;
  881.  
  882. UnityFractal()
  883. {
  884. #ifndef XFRACT
  885.    /* brought to you by Mark Peterson - you won't find this in any fractal
  886.       books unless they saw it here first - Mark invented it! */
  887.    XXOne = multiply(lold.x, lold.x, bitshift) + multiply(lold.y, lold.y, bitshift);
  888.    if((XXOne > FgTwo) || (labs(XXOne - FgOne) < delmin))
  889.       return(1);
  890.    lold.y = multiply(FgTwo - XXOne, lold.x, bitshift);
  891.    lold.x = multiply(FgTwo - XXOne, lold.y, bitshift);
  892.    lnew=lold;  /* TW added this line */
  893.    return(0);
  894. #endif
  895. }
  896.  
  897. #define XXOne new.x
  898.  
  899. UnityfpFractal()
  900. {
  901.    /* brought to you by Mark Peterson - you won't find this in any fractal
  902.       books unless they saw it here first - Mark invented it! */
  903.  
  904.    XXOne = sqr(old.x) + sqr(old.y);
  905.    if((XXOne > 2.0) || (fabs(XXOne - 1.0) < ddelmin))
  906.       return(1);
  907.    old.y = (2.0 - XXOne)* old.x;
  908.    old.x = (2.0 - XXOne)* old.y;
  909.    new=old;  /* TW added this line */
  910.    return(0);
  911. }
  912.  
  913. #undef XXOne
  914.  
  915. Mandel4Fractal()
  916. {
  917.    /* By writing this code, Bert has left behind the excuse "don't
  918.       know what a fractal is, just know how to make'em go fast".
  919.       Bert is hereby declared a bonafide fractal expert! Supposedly
  920.       this routine calculates the Mandelbrot/Julia set based on the
  921.       polynomial z**4 + lambda, but I wouldn't know -- can't follow
  922.       all that integer math speedup stuff - Tim */
  923.  
  924.    /* first, compute (x + iy)**2 */
  925.    lnew.x  = ltempsqrx - ltempsqry;
  926.    lnew.y = multiply(lold.x, lold.y, bitshiftless1);
  927.    if (longbailout()) return(1);
  928.  
  929.    /* then, compute ((x + iy)**2)**2 + lambda */
  930.    lnew.x  = ltempsqrx - ltempsqry + longparm->x;
  931.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  932.    return(longbailout());
  933. }
  934.  
  935. floatZtozPluszpwrFractal()
  936. {
  937.    cpower(&old,(int)param[2],&new);
  938.    old = ComplexPower(old,old);
  939.    new.x = new.x + old.x +floatparm->x;
  940.    new.y = new.y + old.y +floatparm->y;
  941.    return(floatbailout());
  942. }
  943.  
  944. longZpowerFractal()
  945. {
  946. #ifndef XFRACT
  947.    if(lcpower(&lold,c_exp,&lnew,bitshift))
  948.       lnew.x = lnew.y = 8L<<bitshift;
  949.    lnew.x += longparm->x;
  950.    lnew.y += longparm->y;
  951.    return(longbailout());
  952. #endif
  953. }
  954.  
  955. longCmplxZpowerFractal()
  956. {
  957. #ifndef XFRACT
  958.    struct complex x, y;
  959.  
  960.    x.x = (double)lold.x / fudge;
  961.    x.y = (double)lold.y / fudge;
  962.    y.x = (double)lparm2.x / fudge;
  963.    y.y = (double)lparm2.y / fudge;
  964.    x = ComplexPower(x, y);
  965.    if(fabs(x.x) < fgLimit && fabs(x.y) < fgLimit) {
  966.       lnew.x = (long)(x.x * fudge);
  967.       lnew.y = (long)(x.y * fudge);
  968.    }
  969.    else
  970.       overflow = 1;
  971.    lnew.x += longparm->x;
  972.    lnew.y += longparm->y;
  973.    return(longbailout());
  974. #endif
  975. }
  976.  
  977. floatZpowerFractal()
  978. {
  979.    cpower(&old,c_exp,&new);
  980.    new.x += floatparm->x;
  981.    new.y += floatparm->y;
  982.    return(floatbailout());
  983. }
  984.  
  985. floatCmplxZpowerFractal()
  986. {
  987.    new = ComplexPower(old, parm2);
  988.    new.x += floatparm->x;
  989.    new.y += floatparm->y;
  990.    return(floatbailout());
  991. }
  992.  
  993. Barnsley3Fractal()
  994. {
  995.    /* An unnamed Mandelbrot/Julia function from "Fractals
  996.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  997.  
  998.    /* calculate intermediate products */
  999.    oldxinitx   = multiply(lold.x, lold.x, bitshift);
  1000.    oldyinity   = multiply(lold.y, lold.y, bitshift);
  1001.    oldxinity   = multiply(lold.x, lold.y, bitshift);
  1002.  
  1003.    /* orbit calculation */
  1004.    if(lold.x > 0)
  1005.    {
  1006.       lnew.x = oldxinitx   - oldyinity - fudge;
  1007.       lnew.y = oldxinity << 1;
  1008.    }
  1009.    else
  1010.    {
  1011.       lnew.x = oldxinitx - oldyinity - fudge
  1012.        + multiply(longparm->x,lold.x,bitshift);
  1013.       lnew.y = oldxinity <<1;
  1014.  
  1015.       /* This term added by Tim Wegner to make dependent on the
  1016.      imaginary part of the parameter. (Otherwise Mandelbrot
  1017.      is uninteresting. */
  1018.       lnew.y += multiply(longparm->y,lold.x,bitshift);
  1019.    }
  1020.    return(longbailout());
  1021. }
  1022.  
  1023. Barnsley3FPFractal()
  1024. {
  1025.    /* An unnamed Mandelbrot/Julia function from "Fractals
  1026.    Everywhere" by Michael Barnsley, p. 292, example 4.1 */
  1027.  
  1028.  
  1029.    /* calculate intermediate products */
  1030.    foldxinitx  = old.x * old.x;
  1031.    foldyinity  = old.y * old.y;
  1032.    foldxinity  = old.x * old.y;
  1033.  
  1034.    /* orbit calculation */
  1035.    if(old.x > 0)
  1036.    {
  1037.       new.x = foldxinitx - foldyinity - 1.0;
  1038.       new.y = foldxinity * 2;
  1039.    }
  1040.    else
  1041.    {
  1042.       new.x = foldxinitx - foldyinity -1.0 + floatparm->x * old.x;
  1043.       new.y = foldxinity * 2;
  1044.  
  1045.       /* This term added by Tim Wegner to make dependent on the
  1046.      imaginary part of the parameter. (Otherwise Mandelbrot
  1047.      is uninteresting. */
  1048.       new.y += floatparm->y * old.x;
  1049.    }
  1050.    return(floatbailout());
  1051. }
  1052.  
  1053. TrigPlusZsquaredFractal()
  1054. {
  1055. #ifndef XFRACT
  1056.    /* From Scientific American, July 1989 */
  1057.    /* A Biomorph              */
  1058.    /* z(n+1) = trig(z(n))+z(n)**2+C      */
  1059.    LCMPLXtrig0(lold,lnew);
  1060.    lnew.x += ltempsqrx - ltempsqry + longparm->x;
  1061.    lnew.y += multiply(lold.x, lold.y, bitshiftless1) + longparm->y;
  1062.    return(longbailout());
  1063. #endif
  1064. }
  1065.  
  1066. TrigPlusZsquaredfpFractal()
  1067. {
  1068.    /* From Scientific American, July 1989 */
  1069.    /* A Biomorph              */
  1070.    /* z(n+1) = trig(z(n))+z(n)**2+C      */
  1071.  
  1072.    CMPLXtrig0(old,new);
  1073.    new.x += tempsqrx - tempsqry + floatparm->x;
  1074.    new.y += 2.0 * old.x * old.y + floatparm->y;
  1075.    return(floatbailout());
  1076. }
  1077.  
  1078. Richard8fpFractal()
  1079. {
  1080.    /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
  1081.    CMPLXtrig0(old,new);
  1082. /*   CMPLXtrig1(*floatparm,tmp); */
  1083.    new.x += tmp.x;
  1084.    new.y += tmp.y;
  1085.    return(floatbailout());
  1086. }
  1087.  
  1088. Richard8Fractal()
  1089. {
  1090. #ifndef XFRACT
  1091.    /*  Richard8 {c = z = pixel: z=sin(z)+sin(pixel),|z|<=50} */
  1092.    LCMPLXtrig0(lold,lnew);
  1093. /*   LCMPLXtrig1(*longparm,ltmp); */
  1094.    lnew.x += ltmp.x;
  1095.    lnew.y += ltmp.y;
  1096.    return(longbailout());
  1097. #endif
  1098. }
  1099.  
  1100. PopcornFractal()
  1101. {
  1102.    extern int row;
  1103.    tmp = old;
  1104.    tmp.x *= 3.0;
  1105.    tmp.y *= 3.0;
  1106.    FPUsincos(&tmp.x,&sinx,&cosx);
  1107.    FPUsincos(&tmp.y,&siny,&cosy);
  1108.    tmp.x = sinx/cosx + old.x;
  1109.    tmp.y = siny/cosy + old.y;
  1110.    FPUsincos(&tmp.x,&sinx,&cosx);
  1111.    FPUsincos(&tmp.y,&siny,&cosy);
  1112.    new.x = old.x - parm.x*siny;
  1113.    new.y = old.y - parm.x*sinx;
  1114.    /*
  1115.    new.x = old.x - parm.x*sin(old.y+tan(3*old.y));
  1116.    new.y = old.y - parm.x*sin(old.x+tan(3*old.x));
  1117.    */
  1118.    if(plot == noplot)
  1119.    {
  1120.       plot_orbit(new.x,new.y,1+row%colors);
  1121.       old = new;
  1122.    }
  1123.    else
  1124.    /* FLOATBAILOUT(); */
  1125.    /* PB The above line was weird, not what it seems to be!  But, bracketing
  1126.      it or always doing it (either of which seem more likely to be what
  1127.      was intended) changes the image for the worse, so I'm not touching it.
  1128.      Same applies to int form in next routine. */
  1129.    /* PB later: recoded inline, still leaving it weird */
  1130.       tempsqrx = sqr(new.x);
  1131.    tempsqry = sqr(new.y);
  1132.    if((magnitude = tempsqrx + tempsqry) >= rqlim) return(1);
  1133.    old = new;
  1134.    return(0);
  1135. }
  1136.  
  1137. LPopcornFractal()
  1138. {
  1139. #ifndef XFRACT
  1140.    extern int row;
  1141.    ltmp = lold;
  1142.    ltmp.x *= 3L;
  1143.    ltmp.y *= 3L;
  1144.    LTRIGARG(ltmp.x);
  1145.    LTRIGARG(ltmp.y);
  1146.    SinCos086(ltmp.x,&lsinx,&lcosx);
  1147.    SinCos086(ltmp.y,&lsiny,&lcosy);
  1148.    ltmp.x = divide(lsinx,lcosx,bitshift) + lold.x;
  1149.    ltmp.y = divide(lsiny,lcosy,bitshift) + lold.y;
  1150.    LTRIGARG(ltmp.x);
  1151.    LTRIGARG(ltmp.y);
  1152.    SinCos086(ltmp.x,&lsinx,&lcosx);
  1153.    SinCos086(ltmp.y,&lsiny,&lcosy);
  1154.    lnew.x = lold.x - multiply(lparm.x,lsiny,bitshift);
  1155.    lnew.y = lold.y - multiply(lparm.x,lsinx,bitshift);
  1156.    if(plot == noplot)
  1157.    {
  1158.       iplot_orbit(lnew.x,lnew.y,1+row%colors);
  1159.       lold = lnew;
  1160.    }
  1161.    else
  1162.       LONGBAILOUT();
  1163.    /* PB above still the old way, is weird, see notes in FP popcorn case */
  1164.    return(0);
  1165. #endif
  1166. }
  1167.  
  1168. int MarksCplxMand(void)
  1169. {
  1170.    tmp.x = tempsqrx - tempsqry;
  1171.    tmp.y = 2*old.x*old.y;
  1172.    FPUcplxmul(&tmp, &Coefficient, &new);
  1173.    new.x += floatparm->x;
  1174.    new.y += floatparm->y;
  1175.    return(floatbailout());
  1176. }
  1177.  
  1178. int SpiderfpFractal(void)
  1179. {
  1180.    /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
  1181.    new.x = tempsqrx - tempsqry + tmp.x;
  1182.    new.y = 2 * old.x * old.y + tmp.y;
  1183.    tmp.x = tmp.x/2 + new.x;
  1184.    tmp.y = tmp.y/2 + new.y;
  1185.    return(floatbailout());
  1186. }
  1187.  
  1188. SpiderFractal(void)
  1189. {
  1190. #ifndef XFRACT
  1191.    /* Spider(XAXIS) { c=z=pixel: z=z*z+c; c=c/2+z, |z|<=4 } */
  1192.    lnew.x  = ltempsqrx - ltempsqry + ltmp.x;
  1193.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y;
  1194.    ltmp.x = (ltmp.x >> 1) + lnew.x;
  1195.    ltmp.y = (ltmp.y >> 1) + lnew.y;
  1196.    return(longbailout());
  1197. #endif
  1198. }
  1199.  
  1200. TetratefpFractal()
  1201. {
  1202.    /* Tetrate(XAXIS) { c=z=pixel: z=c^z, |z|<=(P1+3) } */
  1203.    new = ComplexPower(*floatparm,old);
  1204.    return(floatbailout());
  1205. }
  1206.  
  1207. ZXTrigPlusZFractal()
  1208. {
  1209. #ifndef XFRACT
  1210.    /* z = (p1*z*trig(z))+p2*z */
  1211.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)         */
  1212.    LCMPLXmult(lparm,ltmp,ltmp);      /* ltmp  = p1*trig(old)          */
  1213.    LCMPLXmult(lold,ltmp,ltmp2);      /* ltmp2 = p1*old*trig(old)      */
  1214.    LCMPLXmult(lparm2,lold,ltmp);     /* ltmp  = p2*old              */
  1215.    LCMPLXadd(ltmp2,ltmp,lnew);         /* lnew  = p1*trig(old) + p2*old */
  1216.    return(longbailout());
  1217. #endif
  1218. }
  1219.  
  1220. ScottZXTrigPlusZFractal()
  1221. {
  1222. #ifndef XFRACT
  1223.    /* z = (z*trig(z))+z */
  1224.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)       */
  1225.    LCMPLXmult(lold,ltmp,lnew);         /* lnew  = old*trig(old)    */
  1226.    LCMPLXadd(lnew,lold,lnew);         /* lnew  = trig(old) + old */
  1227.    return(longbailout());
  1228. #endif
  1229. }
  1230.  
  1231. SkinnerZXTrigSubZFractal()
  1232. {
  1233. #ifndef XFRACT
  1234.    /* z = (z*trig(z))-z */
  1235.    LCMPLXtrig0(lold,ltmp);        /* ltmp  = trig(old)       */
  1236.    LCMPLXmult(lold,ltmp,lnew);         /* lnew  = old*trig(old)    */
  1237.    LCMPLXsub(lnew,lold,lnew);         /* lnew  = trig(old) - old */
  1238.    return(longbailout());
  1239. #endif
  1240. }
  1241.  
  1242. ZXTrigPlusZfpFractal()
  1243. {
  1244.    /* z = (p1*z*trig(z))+p2*z */
  1245.    CMPLXtrig0(old,tmp);      /* tmp  = trig(old)         */
  1246.    CMPLXmult(parm,tmp,tmp);     /* tmp  = p1*trig(old)      */
  1247.    CMPLXmult(old,tmp,tmp2);     /* tmp2 = p1*old*trig(old)     */
  1248.    CMPLXmult(parm2,old,tmp);     /* tmp  = p2*old         */
  1249.    CMPLXadd(tmp2,tmp,new);     /* new  = p1*trig(old) + p2*old */
  1250.    return(floatbailout());
  1251. }
  1252.  
  1253. ScottZXTrigPlusZfpFractal()
  1254. {
  1255.    /* z = (z*trig(z))+z */
  1256.    CMPLXtrig0(old,tmp);     /* tmp    = trig(old)      */
  1257.    CMPLXmult(old,tmp,new);     /* new  = old*trig(old)   */
  1258.    CMPLXadd(new,old,new);     /* new  = trig(old) + old */
  1259.    return(floatbailout());
  1260. }
  1261.  
  1262. SkinnerZXTrigSubZfpFractal()
  1263. {
  1264.    /* z = (z*trig(z))-z */
  1265.    CMPLXtrig0(old,tmp);     /* tmp    = trig(old)      */
  1266.    CMPLXmult(old,tmp,new);     /* new  = old*trig(old)   */
  1267.    CMPLXsub(new,old,new);     /* new  = trig(old) - old */
  1268.    return(floatbailout());
  1269. }
  1270.  
  1271. Sqr1overTrigFractal()
  1272. {
  1273. #ifndef XFRACT
  1274.    /* z = sqr(1/trig(z)) */
  1275.    LCMPLXtrig0(lold,lold);
  1276.    LCMPLXrecip(lold,lold);
  1277.    LCMPLXsqr(lold,lnew);
  1278.    return(longbailout());
  1279. #endif
  1280. }
  1281.  
  1282. Sqr1overTrigfpFractal()
  1283. {
  1284.    /* z = sqr(1/trig(z)) */
  1285.    CMPLXtrig0(old,old);
  1286.    CMPLXrecip(old,old);
  1287.    CMPLXsqr(old,new);
  1288.    return(floatbailout());
  1289. }
  1290.  
  1291. TrigPlusTrigFractal()
  1292. {
  1293. #ifndef XFRACT
  1294.    /* z = trig(0,z)*p1+trig1(z)*p2 */
  1295.    LCMPLXtrig0(lold,ltmp);
  1296.    LCMPLXmult(lparm,ltmp,ltmp);
  1297.    LCMPLXtrig1(lold,ltmp2);
  1298.    LCMPLXmult(lparm2,ltmp2,lold);
  1299.    LCMPLXadd(ltmp,lold,lnew);
  1300.    return(longbailout());
  1301. #endif
  1302. }
  1303.  
  1304. TrigPlusTrigfpFractal()
  1305. {
  1306.    /* z = trig0(z)*p1+trig1(z)*p2 */
  1307.    CMPLXtrig0(old,tmp);
  1308.    CMPLXmult(parm,tmp,tmp);
  1309.    CMPLXtrig1(old,old);
  1310.    CMPLXmult(parm2,old,old);
  1311.    CMPLXadd(tmp,old,new);
  1312.    return(floatbailout());
  1313. }
  1314.  
  1315. /* The following four fractals are based on the idea of parallel
  1316.    or alternate calculations.  The shift is made when the mod
  1317.    reaches a given value.  JCO  5/6/92 */
  1318.  
  1319. LambdaTrigOrTrigFractal()
  1320. {
  1321. #ifndef XFRACT
  1322.    /* z = trig0(z)*p1 if mod(old) < p2.x and
  1323.           trig1(z)*p1 if mod(old) >= p2.x */
  1324.    if ((LCMPLXmod(lold)) < lparm2.x){
  1325.      LCMPLXtrig0(lold,ltmp);
  1326.      LCMPLXmult(lparm,ltmp,lnew);}
  1327.    else{
  1328.      LCMPLXtrig1(lold,ltmp);
  1329.      LCMPLXmult(lparm,ltmp,lnew);}
  1330.    return(longbailout());
  1331. #endif
  1332. }
  1333.  
  1334. LambdaTrigOrTrigfpFractal()
  1335. {
  1336.    /* z = trig0(z)*p1 if mod(old) < p2.x and
  1337.           trig1(z)*p1 if mod(old) >= p2.x */
  1338.    if (CMPLXmod(old) < parm2.x){
  1339.      CMPLXtrig0(old,old);
  1340.      CMPLXmult(parm,old,new);}
  1341.    else{
  1342.      CMPLXtrig1(old,old);
  1343.      CMPLXmult(parm,old,new);}
  1344.    return(floatbailout());
  1345. }
  1346.  
  1347. JuliaTrigOrTrigFractal()
  1348. {
  1349. #ifndef XFRACT
  1350.    /* z = trig0(z)+p1 if mod(old) < p2.x and
  1351.           trig1(z)+p1 if mod(old) >= p2.x */
  1352.    if (LCMPLXmod(lold) < lparm2.x){
  1353.      LCMPLXtrig0(lold,ltmp);
  1354.      LCMPLXadd(lparm,ltmp,lnew);}
  1355.    else{
  1356.      LCMPLXtrig1(lold,ltmp);
  1357.      LCMPLXadd(lparm,ltmp,lnew);}
  1358.    return(longbailout());
  1359. #endif
  1360. }
  1361.  
  1362. JuliaTrigOrTrigfpFractal()
  1363. {
  1364.    /* z = trig0(z)+p1 if mod(old) < p2.x and
  1365.           trig1(z)+p1 if mod(old) >= p2.x */
  1366.    if (CMPLXmod(old) < parm2.x){
  1367.      CMPLXtrig0(old,old);
  1368.      CMPLXadd(parm,old,new);}
  1369.    else{
  1370.      CMPLXtrig1(old,old);
  1371.      CMPLXadd(parm,old,new);}
  1372.    return(floatbailout());
  1373. }
  1374.  
  1375. ManlamTrigOrTrigFractal()
  1376. {  /* Psuedo Mandelbrot set corresponding to LambdaTrigOrTrigFractal */
  1377. #ifndef XFRACT
  1378.    /* z = trig0(z)*pixel if mod(old) < p2.x and
  1379.           trig1(z)*pixel if mod(old) >= p2.x */
  1380.    if (LCMPLXmod(lold) < lparm2.x){
  1381.      LCMPLXtrig0(lold,ltmp);
  1382.      LCMPLXmult(linit,ltmp,lnew);}
  1383.    else{
  1384.      LCMPLXtrig1(lold,ltmp);
  1385.      LCMPLXmult(linit,ltmp,lnew);}
  1386.    return(longbailout());
  1387. #endif
  1388. }
  1389.  
  1390. ManlamTrigOrTrigfpFractal()
  1391. {  /* Psuedo Mandelbrot set corresponding to LambdaTrigOrTrigfpFractal */
  1392.    /* z = trig0(z)*pixel if mod(old) < p2.x and
  1393.           trig1(z)*pixel if mod(old) >= p2.x */
  1394.    if (CMPLXmod(old) < parm2.x){
  1395.      CMPLXtrig0(old,tmp);
  1396.      CMPLXmult(init,tmp,new);}
  1397.    else{
  1398.      CMPLXtrig1(old,tmp);
  1399.      CMPLXmult(init,tmp,new);}
  1400.    return(floatbailout());
  1401. }
  1402.  
  1403. MandelTrigOrTrigFractal()
  1404. {
  1405. #ifndef XFRACT
  1406.    /* z = trig0(z)+pixel if mod(old) < p2.x and
  1407.           trig1(z)+pixel if mod(old) >= p2.x */
  1408.    if (LCMPLXmod(lold) < lparm2.x){
  1409.      LCMPLXtrig0(lold,ltmp);
  1410.      LCMPLXadd(linit,ltmp,lnew);}
  1411.    else{
  1412.      LCMPLXtrig1(lold,ltmp);
  1413.      LCMPLXadd(linit,ltmp,lnew);}
  1414.    return(longbailout());
  1415. #endif
  1416. }
  1417.  
  1418. MandelTrigOrTrigfpFractal()
  1419. {
  1420.    /* z = trig0(z)+pixel if mod(old) < p2.x and
  1421.           trig1(z)+pixel if mod(old) >= p2.x */
  1422.    if (CMPLXmod(old) < parm2.x){
  1423.      CMPLXtrig0(old,tmp);
  1424.      CMPLXadd(init,tmp,new);}
  1425.    else{
  1426.      CMPLXtrig1(old,tmp);
  1427.      CMPLXadd(init,tmp,new);}
  1428.    return(floatbailout());
  1429. }
  1430.  
  1431. static int AplusOne, Ap1deg;
  1432. static struct MP mpAplusOne, mpAp1deg, mpdegree, mptmpparmy;
  1433.  
  1434. int MPCHalleyFractal()
  1435. {
  1436. #ifndef XFRACT
  1437.    /*  X(X^a - 1) = 0, Halley Map */
  1438.    /*  a = parm.x,  relaxation coeff. = parm.y,  epsilon = parm2.x  */
  1439.  
  1440. int ihal;
  1441. struct MPC mpcXtoAlessOne, mpcXtoA;
  1442. struct MPC mpcXtoAplusOne; /* a-1, a, a+1 */
  1443. struct MPC mpcFX, mpcF1prime, mpcF2prime, mpcHalnumer1;
  1444. struct MPC mpcHalnumer2, mpcHaldenom, mpctmp;
  1445.  
  1446.    mpcXtoAlessOne = mpcold;
  1447.    for(ihal=2; ihal<degree; ihal++) {
  1448.      mpcXtoAlessOne.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
  1449.      mpcXtoAlessOne.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
  1450.    }
  1451.    mpcXtoA.x = *pMPsub(*pMPmul(mpcXtoAlessOne.x,mpcold.x),*pMPmul(mpcXtoAlessOne.y,mpcold.y));
  1452.    mpcXtoA.y = *pMPadd(*pMPmul(mpcXtoAlessOne.x,mpcold.y),*pMPmul(mpcXtoAlessOne.y,mpcold.x));
  1453.    mpcXtoAplusOne.x = *pMPsub(*pMPmul(mpcXtoA.x,mpcold.x),*pMPmul(mpcXtoA.y,mpcold.y));
  1454.    mpcXtoAplusOne.y = *pMPadd(*pMPmul(mpcXtoA.x,mpcold.y),*pMPmul(mpcXtoA.y,mpcold.x));
  1455.  
  1456.    mpcFX.x = *pMPsub(mpcXtoAplusOne.x, mpcold.x);
  1457.    mpcFX.y = *pMPsub(mpcXtoAplusOne.y, mpcold.y); /* FX = X^(a+1) - X  = F */
  1458.  
  1459.    mpcF2prime.x = *pMPmul(mpAp1deg, mpcXtoAlessOne.x); /* mpAp1deg in setup */
  1460.    mpcF2prime.y = *pMPmul(mpAp1deg, mpcXtoAlessOne.y);        /* F" */
  1461.  
  1462.    mpcF1prime.x = *pMPsub(*pMPmul(mpAplusOne, mpcXtoA.x), mpone);
  1463.    mpcF1prime.y = *pMPmul(mpAplusOne, mpcXtoA.y);                   /*  F'  */
  1464.  
  1465.    mpcHalnumer1.x = *pMPsub(*pMPmul(mpcF2prime.x,mpcFX.x),*pMPmul(mpcF2prime.y,mpcFX.y));
  1466.    mpcHalnumer1.y = *pMPadd(*pMPmul(mpcF2prime.x,mpcFX.y),*pMPmul(mpcF2prime.y,mpcFX.x));
  1467.    /*  F * F"  */
  1468.  
  1469.    mpcHaldenom.x = *pMPadd(mpcF1prime.x, mpcF1prime.x);
  1470.    mpcHaldenom.y = *pMPadd(mpcF1prime.y, mpcF1prime.y);      /*  2 * F'  */
  1471.  
  1472.    mpcHalnumer1 = MPCdiv(mpcHalnumer1, mpcHaldenom);        /*  F"F/2F'  */
  1473.    mpcHalnumer2.x = *pMPsub(mpcF1prime.x, mpcHalnumer1.x);
  1474.    mpcHalnumer2.y = *pMPsub(mpcF1prime.y, mpcHalnumer1.y); /*  F' - F"F/2F'  */
  1475.    mpcHalnumer2 = MPCdiv(mpcFX, mpcHalnumer2);
  1476.  
  1477.    mpctmp.x = *pMPmul(mptmpparmy,mpcHalnumer2.x); /* mptmpparmy is */
  1478.    mpctmp.y = *pMPmul(mptmpparmy,mpcHalnumer2.y); /* relaxation coef. */
  1479.    mpcnew.x = *pMPsub(mpcold.x, mpctmp.x);
  1480.    mpcnew.y = *pMPsub(mpcold.y, mpctmp.y);
  1481.  
  1482.    return(MPCHalleybailout());
  1483. #endif
  1484. }
  1485.  
  1486. HalleyFractal()
  1487. {
  1488.    /*  X(X^a - 1) = 0, Halley Map */
  1489.    /*  a = parm.x = degree, relaxation coeff. = parm.y, epsilon = parm2.x  */
  1490.  
  1491. int ihal;
  1492. CMPLX XtoAlessOne, XtoA, XtoAplusOne; /* a-1, a, a+1 */
  1493. CMPLX FX, F1prime, F2prime, Halnumer1, Halnumer2, Haldenom;
  1494.  
  1495.    XtoAlessOne = old;
  1496.    for(ihal=2; ihal<degree; ihal++) {
  1497.      CMPLXmult(old, XtoAlessOne, XtoAlessOne);
  1498.    }
  1499.    CMPLXmult(old, XtoAlessOne, XtoA);
  1500.    CMPLXmult(old, XtoA, XtoAplusOne);
  1501.  
  1502.    CMPLXsub(XtoAplusOne, old, FX);        /* FX = X^(a+1) - X  = F */
  1503.    F2prime.x = Ap1deg * XtoAlessOne.x; /* Ap1deg in setup */
  1504.    F2prime.y = Ap1deg * XtoAlessOne.y;        /* F" */
  1505.  
  1506.    F1prime.x = AplusOne * XtoA.x - 1.0;
  1507.    F1prime.y = AplusOne * XtoA.y;                             /*  F'  */
  1508.  
  1509.    CMPLXmult(F2prime, FX, Halnumer1);                      /*  F * F"  */
  1510.    Haldenom.x = F1prime.x + F1prime.x;
  1511.    Haldenom.y = F1prime.y + F1prime.y;                     /*  2 * F'  */
  1512.  
  1513.    complex_div(Halnumer1, Haldenom, &Halnumer1);          /*  F"F/2F'  */
  1514.    CMPLXsub(F1prime, Halnumer1, Halnumer2);          /*  F' - F"F/2F'  */
  1515.    complex_div(FX, Halnumer2, &Halnumer2);
  1516.    new.x = old.x - (parm.y * Halnumer2.x); /* parm.y is relaxation coef. */
  1517.    new.y = old.y - (parm.y * Halnumer2.y);
  1518.  
  1519.    return(Halleybailout());
  1520. }
  1521.  
  1522. ScottTrigPlusTrigFractal()
  1523. {
  1524. #ifndef XFRACT
  1525.    /* z = trig0(z)+trig1(z) */
  1526.    LCMPLXtrig0(lold,ltmp);
  1527.    LCMPLXtrig1(lold,lold);
  1528.    LCMPLXadd(ltmp,lold,lnew);
  1529.    return(longbailout());
  1530. #endif
  1531. }
  1532.  
  1533. ScottTrigPlusTrigfpFractal()
  1534. {
  1535.    /* z = trig0(z)+trig1(z) */
  1536.    CMPLXtrig0(old,tmp);
  1537.    CMPLXtrig1(old,tmp2);
  1538.    CMPLXadd(tmp,tmp2,new);
  1539.    return(floatbailout());
  1540. }
  1541.  
  1542. SkinnerTrigSubTrigFractal()
  1543. {
  1544. #ifndef XFRACT
  1545.    /* z = trig(0,z)-trig1(z) */
  1546.    LCMPLXtrig0(lold,ltmp);
  1547.    LCMPLXtrig1(lold,ltmp2);
  1548.    LCMPLXsub(ltmp,ltmp2,lnew);
  1549.    return(longbailout());
  1550. #endif
  1551. }
  1552.  
  1553. SkinnerTrigSubTrigfpFractal()
  1554. {
  1555.    /* z = trig0(z)-trig1(z) */
  1556.    CMPLXtrig0(old,tmp);
  1557.    CMPLXtrig1(old,tmp2);
  1558.    CMPLXsub(tmp,tmp2,new);
  1559.    return(floatbailout());
  1560. }
  1561.  
  1562.  
  1563. TrigXTrigfpFractal()
  1564. {
  1565.    /* z = trig0(z)*trig1(z) */
  1566.    CMPLXtrig0(old,tmp);
  1567.    CMPLXtrig1(old,old);
  1568.    CMPLXmult(tmp,old,new);
  1569.    return(floatbailout());
  1570. }
  1571.  
  1572. TrigXTrigFractal()
  1573. {
  1574. #ifndef XFRACT
  1575.    LCMPLX ltmp2;
  1576.    /* z = trig0(z)*trig1(z) */
  1577.    LCMPLXtrig0(lold,ltmp);
  1578.    LCMPLXtrig1(lold,ltmp2);
  1579.    LCMPLXmult(ltmp,ltmp2,lnew);
  1580.    if(overflow)
  1581.       TryFloatFractal(TrigXTrigfpFractal);
  1582.    return(longbailout());
  1583. #endif
  1584. }
  1585.  
  1586.  /* call float version of fractal if integer math overflow */
  1587. TryFloatFractal(int (*fpFractal)())
  1588. {
  1589.    overflow=0;
  1590.    /* lold had better not be changed! */
  1591.    old.x = lold.x; old.x /= fudge;
  1592.    old.y = lold.y; old.y /= fudge;
  1593.    tempsqrx = sqr(old.x);
  1594.    tempsqry = sqr(old.y);
  1595.    fpFractal();
  1596.    lnew.x = new.x/fudge;
  1597.    lnew.y = new.y/fudge;
  1598.    return(0);
  1599. }
  1600.  
  1601. /********************************************************************/
  1602. /*  Next six orbit functions are one type - extra functions are     */
  1603. /*    special cases written for speed.                    */
  1604. /********************************************************************/
  1605.  
  1606. TrigPlusSqrFractal() /* generalization of Scott and Skinner types */
  1607. {
  1608. #ifndef XFRACT
  1609.    /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
  1610.    LCMPLXtrig0(lold,ltmp);     /* ltmp = trig(lold)               */
  1611.    LCMPLXmult(lparm,ltmp,lnew); /* lnew = lparm*trig(lold)            */
  1612.    LCMPLXsqr_old(ltmp);     /* ltmp = sqr(lold)                */
  1613.    LCMPLXmult(lparm2,ltmp,ltmp);/* ltmp = lparm2*sqr(lold)            */
  1614.    LCMPLXadd(lnew,ltmp,lnew);    /* lnew = lparm*trig(lold)+lparm2*sqr(lold) */
  1615.    return(longbailout());
  1616. #endif
  1617. }
  1618.  
  1619. TrigPlusSqrfpFractal() /* generalization of Scott and Skinner types */
  1620. {
  1621.    /* { z=pixel: z=(p1,p2)*trig(z)+(p3,p4)*sqr(z), |z|<BAILOUT } */
  1622.    CMPLXtrig0(old,tmp);     /* tmp = trig(old)               */
  1623.    CMPLXmult(parm,tmp,new); /* new = parm*trig(old)           */
  1624.    CMPLXsqr_old(tmp);         /* tmp = sqr(old)                */
  1625.    CMPLXmult(parm2,tmp,tmp2); /* tmp = parm2*sqr(old)             */
  1626.    CMPLXadd(new,tmp2,new);    /* new = parm*trig(old)+parm2*sqr(old) */
  1627.    return(floatbailout());
  1628. }
  1629.  
  1630. ScottTrigPlusSqrFractal()
  1631. {
  1632. #ifndef XFRACT
  1633.    /*  { z=pixel: z=trig(z)+sqr(z), |z|<BAILOUT } */
  1634.    LCMPLXtrig0(lold,lnew);    /* lnew = trig(lold)         */
  1635.    LCMPLXsqr_old(ltmp);        /* lold = sqr(lold)          */
  1636.    LCMPLXadd(ltmp,lnew,lnew);  /* lnew = trig(lold)+sqr(lold) */
  1637.    return(longbailout());
  1638. #endif
  1639. }
  1640.  
  1641. ScottTrigPlusSqrfpFractal() /* float version */
  1642. {
  1643.    /* { z=pixel: z=sin(z)+sqr(z), |z|<BAILOUT } */
  1644.    CMPLXtrig0(old,new);       /* new = trig(old)      */
  1645.    CMPLXsqr_old(tmp);           /* tmp = sqr(old)       */
  1646.    CMPLXadd(new,tmp,new);      /* new = trig(old)+sqr(old) */
  1647.    return(floatbailout());
  1648. }
  1649.  
  1650. SkinnerTrigSubSqrFractal()
  1651. {
  1652. #ifndef XFRACT
  1653.    /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT }           */
  1654.    LCMPLXtrig0(lold,lnew);    /* lnew = trig(lold)         */
  1655.    LCMPLXsqr_old(ltmp);        /* lold = sqr(lold)          */
  1656.    LCMPLXsub(lnew,ltmp,lnew);  /* lnew = trig(lold)-sqr(lold) */
  1657.    return(longbailout());
  1658. #endif
  1659. }
  1660.  
  1661. SkinnerTrigSubSqrfpFractal()
  1662. {
  1663.    /* { z=pixel: z=sin(z)-sqr(z), |z|<BAILOUT } */
  1664.    CMPLXtrig0(old,new);       /* new = trig(old) */
  1665.    CMPLXsqr_old(tmp);           /* old = sqr(old)  */
  1666.    CMPLXsub(new,tmp,new);      /* new = trig(old)-sqr(old) */
  1667.    return(floatbailout());
  1668. }
  1669.  
  1670. TrigZsqrdfpFractal()
  1671. {
  1672.    /* { z=pixel: z=trig(z*z), |z|<TEST } */
  1673.    CMPLXsqr_old(tmp);
  1674.    CMPLXtrig0(tmp,new);
  1675.    return(floatbailout());
  1676. }
  1677.  
  1678. TrigZsqrdFractal() /* this doesn't work very well */
  1679. {
  1680. #ifndef XFRACT
  1681.    /* { z=pixel: z=trig(z*z), |z|<TEST } */
  1682.    LCMPLXsqr_old(ltmp);
  1683.    LCMPLXtrig0(ltmp,lnew);
  1684.    if(overflow)
  1685.       TryFloatFractal(TrigZsqrdfpFractal);
  1686.    return(longbailout());
  1687. #endif
  1688. }
  1689.  
  1690. SqrTrigFractal()
  1691. {
  1692. #ifndef XFRACT
  1693.    /* { z=pixel: z=sqr(trig(z)), |z|<TEST} */
  1694.    LCMPLXtrig0(lold,ltmp);
  1695.    LCMPLXsqr(ltmp,lnew);
  1696.    return(longbailout());
  1697. #endif
  1698. }
  1699.  
  1700. SqrTrigfpFractal()
  1701. {
  1702.    /* SZSB(XYAXIS) { z=pixel, TEST=(p1+3): z=sin(z)*sin(z), |z|<TEST} */
  1703.    CMPLXtrig0(old,tmp);
  1704.    CMPLXsqr(tmp,new);
  1705.    return(floatbailout());
  1706. }
  1707.  
  1708.  
  1709. Magnet1Fractal()    /*      Z = ((Z**2 + C - 1)/(2Z + C - 2))**2      */
  1710.   {              /*  In "Beauty of Fractals", code by Kev Allen. */
  1711.     CMPLX top, bot, tmp;
  1712.     double div;
  1713.  
  1714.     top.x = tempsqrx - tempsqry + floatparm->x - 1; /* top = Z**2+C-1 */
  1715.     top.y = old.x * old.y;
  1716.     top.y = top.y + top.y + floatparm->y;
  1717.  
  1718.     bot.x = old.x + old.x + floatparm->x - 2;        /* bot = 2*Z+C-2  */
  1719.     bot.y = old.y + old.y + floatparm->y;
  1720.  
  1721.     div = bot.x*bot.x + bot.y*bot.y;            /* tmp = top/bot  */
  1722.     if (div < FLT_MIN) return(1);
  1723.     tmp.x = (top.x*bot.x + top.y*bot.y)/div;
  1724.     tmp.y = (top.y*bot.x - top.x*bot.y)/div;
  1725.  
  1726.     new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y);        /* Z = tmp**2     */
  1727.     new.y = tmp.x * tmp.y;
  1728.     new.y += new.y;
  1729.  
  1730.     return(floatbailout());
  1731.   }
  1732.  
  1733. Magnet2Fractal()  /* Z = ((Z**3 + 3(C-1)Z + (C-1)(C-2)    ) /     */
  1734.             /*         (3Z**2 + 3(C-2)Z + (C-1)(C-2)+1) )**2  */
  1735.   {            /*     In "Beauty of Fractals", code by Kev Allen.  */
  1736.     CMPLX top, bot, tmp;
  1737.     double div;
  1738.  
  1739.     top.x = old.x * (tempsqrx-tempsqry-tempsqry-tempsqry + T_Cm1.x)
  1740.       - old.y * T_Cm1.y + T_Cm1Cm2.x;
  1741.     top.y = old.y * (tempsqrx+tempsqrx+tempsqrx-tempsqry + T_Cm1.x)
  1742.       + old.x * T_Cm1.y + T_Cm1Cm2.y;
  1743.  
  1744.     bot.x = tempsqrx - tempsqry;
  1745.     bot.x = bot.x + bot.x + bot.x
  1746.       + old.x * T_Cm2.x - old.y * T_Cm2.y
  1747.       + T_Cm1Cm2.x + 1.0;
  1748.     bot.y = old.x * old.y;
  1749.     bot.y += bot.y;
  1750.     bot.y = bot.y + bot.y + bot.y
  1751.       + old.x * T_Cm2.y + old.y * T_Cm2.x
  1752.       + T_Cm1Cm2.y;
  1753.  
  1754.     div = bot.x*bot.x + bot.y*bot.y;            /* tmp = top/bot  */
  1755.     if (div < FLT_MIN) return(1);
  1756.     tmp.x = (top.x*bot.x + top.y*bot.y)/div;
  1757.     tmp.y = (top.y*bot.x - top.x*bot.y)/div;
  1758.  
  1759.     new.x = (tmp.x + tmp.y) * (tmp.x - tmp.y);        /* Z = tmp**2     */
  1760.     new.y = tmp.x * tmp.y;
  1761.     new.y += new.y;
  1762.  
  1763.     return(floatbailout());
  1764.   }
  1765.  
  1766. LambdaTrigFractal()
  1767. {
  1768. #ifndef XFRACT
  1769.    LONGXYTRIGBAILOUT();
  1770.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1771.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1772.    lold = lnew;
  1773.    return(0);
  1774. #endif
  1775. }
  1776.  
  1777. LambdaTrigfpFractal()
  1778. {
  1779.    FLOATXYTRIGBAILOUT();
  1780.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1781.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1782.    old = new;
  1783.    return(0);
  1784. }
  1785.  
  1786. /* bailouts are different for different trig functions */
  1787. LambdaTrigFractal1()
  1788. {
  1789. #ifndef XFRACT
  1790.    LONGTRIGBAILOUT(); /* sin,cos */
  1791.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1792.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1793.    lold = lnew;
  1794.    return(0);
  1795. #endif
  1796. }
  1797.  
  1798. LambdaTrigfpFractal1()
  1799. {
  1800.    FLOATTRIGBAILOUT(); /* sin,cos */
  1801.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1802.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1803.    old = new;
  1804.    return(0);
  1805. }
  1806.  
  1807. LambdaTrigFractal2()
  1808. {
  1809. #ifndef XFRACT
  1810.    LONGHTRIGBAILOUT(); /* sinh,cosh */
  1811.    LCMPLXtrig0(lold,ltmp);         /* ltmp = trig(lold)        */
  1812.    LCMPLXmult(*longparm,ltmp,lnew);   /* lnew = longparm*trig(lold)  */
  1813.    lold = lnew;
  1814.    return(0);
  1815. #endif
  1816. }
  1817.  
  1818. LambdaTrigfpFractal2()
  1819. {
  1820. #ifndef XFRACT
  1821.    FLOATHTRIGBAILOUT(); /* sinh,cosh */
  1822.    CMPLXtrig0(old,tmp);          /* tmp = trig(old)       */
  1823.    CMPLXmult(*floatparm,tmp,new);   /* new = longparm*trig(old)  */
  1824.    old = new;
  1825.    return(0);
  1826. #endif
  1827. }
  1828.  
  1829. ManOWarFractal()
  1830. {
  1831. #ifndef XFRACT
  1832.    /* From Art Matrix via Lee Skinner */
  1833.    lnew.x  = ltempsqrx - ltempsqry + ltmp.x + longparm->x;
  1834.    lnew.y = multiply(lold.x, lold.y, bitshiftless1) + ltmp.y + longparm->y;
  1835.    ltmp = lold;
  1836.    return(longbailout());
  1837. #endif
  1838. }
  1839.  
  1840. ManOWarfpFractal()
  1841. {
  1842.    /* From Art Matrix via Lee Skinner */
  1843.    /* note that fast >= 287 equiv in fracsuba.asm must be kept in step */
  1844.    new.x = tempsqrx - tempsqry + tmp.x + floatparm->x;
  1845.    new.y = 2.0 * old.x * old.y + tmp.y + floatparm->y;
  1846.    tmp = old;
  1847.    return(floatbailout());
  1848. }
  1849.  
  1850. /*
  1851.    MarksMandelPwr (XAXIS) {
  1852.       z = pixel, c = z ^ (z - 1):
  1853.      z = c * sqr(z) + pixel,
  1854.       |z| <= 4
  1855.    }
  1856. */
  1857.  
  1858. MarksMandelPwrfpFractal()
  1859. {
  1860.    CMPLXtrig0(old,new);
  1861.    CMPLXmult(tmp,new,new);
  1862.    new.x += floatparm->x;
  1863.    new.y += floatparm->y;
  1864.    return(floatbailout());
  1865. }
  1866.  
  1867. MarksMandelPwrFractal()
  1868. {
  1869. #ifndef XFRACT
  1870.    LCMPLXtrig0(lold,lnew);
  1871.    LCMPLXmult(ltmp,lnew,lnew);
  1872.    lnew.x += longparm->x;
  1873.    lnew.y += longparm->y;
  1874.    return(longbailout());
  1875. #endif
  1876. }
  1877.  
  1878. /* I was coding Marksmandelpower and failed to use some temporary
  1879.    variables. The result was nice, and since my name is not on any fractal,
  1880.    I thought I would immortalize myself with this error!
  1881.         Tim Wegner */
  1882.  
  1883.  
  1884. TimsErrorfpFractal()
  1885. {
  1886.    CMPLXtrig0(old,new);
  1887.    new.x = new.x * tmp.x - new.y * tmp.y;
  1888.    new.y = new.x * tmp.y - new.y * tmp.x;
  1889.    new.x += floatparm->x;
  1890.    new.y += floatparm->y;
  1891.    return(floatbailout());
  1892. }
  1893. TimsErrorFractal()
  1894. {
  1895. #ifndef XFRACT
  1896.    LCMPLXtrig0(lold,lnew);
  1897.    lnew.x = multiply(lnew.x,ltmp.x,bitshift)-multiply(lnew.y,ltmp.y,bitshift);
  1898.    lnew.y = multiply(lnew.x,ltmp.y,bitshift)-multiply(lnew.y,ltmp.x,bitshift);
  1899.    lnew.x += longparm->x;
  1900.    lnew.y += longparm->y;
  1901.    return(longbailout());
  1902. #endif
  1903. }
  1904.  
  1905. CirclefpFractal()
  1906. {
  1907.    extern int colors;
  1908.    extern int color;
  1909.    int i;
  1910.    i = param[0]*(tempsqrx+tempsqry);
  1911.    color = i&(colors-1);
  1912.    return(1);
  1913. }
  1914. /*
  1915. CirclelongFractal()
  1916. {
  1917.    extern int colors;
  1918.    extern int color;
  1919.    long i;
  1920.    i = multiply(lparm.x,(ltempsqrx+ltempsqry),bitshift);
  1921.    i = i >> bitshift;
  1922.    color = i&(colors-1);
  1923.    return(1);
  1924. }
  1925. */
  1926.  
  1927. /* -------------------------------------------------------------------- */
  1928. /*        Initialization (once per pixel) routines                        */
  1929. /* -------------------------------------------------------------------- */
  1930.  
  1931. #ifdef XFRACT
  1932. /* this code translated to asm - lives in newton.asm */
  1933. /* transform points with reciprocal function */
  1934. void invertz2(CMPLX *z)
  1935. {
  1936.    z->x = dx0[col]+dx1[row];
  1937.    z->y = dy0[row]+dy1[col];
  1938.    z->x -= f_xcenter; z->y -= f_ycenter;  /* Normalize values to center of circle */
  1939.  
  1940.    tempsqrx = sqr(z->x) + sqr(z->y);  /* Get old radius */
  1941.    if(fabs(tempsqrx) > FLT_MIN)
  1942.       tempsqrx = f_radius / tempsqrx;
  1943.    else
  1944.       tempsqrx = FLT_MAX;   /* a big number, but not TOO big */
  1945.    z->x *= tempsqrx;      z->y *= tempsqrx;     /* Perform inversion */
  1946.    z->x += f_xcenter; z->y += f_ycenter; /* Renormalize */
  1947. }
  1948. #endif
  1949.  
  1950. int long_julia_per_pixel()
  1951. {
  1952. #ifndef XFRACT
  1953.    /* integer julia types */
  1954.    /* lambda */
  1955.    /* barnsleyj1 */
  1956.    /* barnsleyj2 */
  1957.    /* sierpinski */
  1958.    if(invert)
  1959.    {
  1960.       /* invert */
  1961.       invertz2(&old);
  1962.  
  1963.       /* watch out for overflow */
  1964.       if(sqr(old.x)+sqr(old.y) >= 127)
  1965.       {
  1966.      old.x = 8;  /* value to bail out in one iteration */
  1967.      old.y = 8;
  1968.       }
  1969.  
  1970.       /* convert to fudged longs */
  1971.       lold.x = old.x*fudge;
  1972.       lold.y = old.y*fudge;
  1973.    }
  1974.    else
  1975.    {
  1976.       lold.x = lx0[col]+lx1[row];
  1977.       lold.y = ly0[row]+ly1[col];
  1978.    }
  1979.    return(0);
  1980. #else
  1981.    printf("Called long_julia_per_pixel\n");
  1982.    exit(0);
  1983. #endif
  1984. }
  1985.  
  1986. int long_richard8_per_pixel()
  1987. {
  1988. #ifndef XFRACT
  1989.     long_mandel_per_pixel();
  1990.     LCMPLXtrig1(*longparm,ltmp);
  1991.     LCMPLXmult(ltmp,lparm2,ltmp);
  1992.     return(1);
  1993. #endif
  1994. }
  1995.  
  1996. int long_mandel_per_pixel()
  1997. {
  1998. #ifndef XFRACT
  1999.    /* integer mandel types */
  2000.    /* barnsleym1 */
  2001.    /* barnsleym2 */
  2002.    linit.x = lx0[col]+lx1[row];
  2003.  
  2004.    if(invert)
  2005.    {
  2006.       /* invert */
  2007.       invertz2(&init);
  2008.  
  2009.       /* watch out for overflow */
  2010.       if(sqr(init.x)+sqr(init.y) >= 127)
  2011.       {
  2012.      init.x = 8;  /* value to bail out in one iteration */
  2013.      init.y = 8;
  2014.       }
  2015.  
  2016.       /* convert to fudged longs */
  2017.       linit.x = init.x*fudge;
  2018.       linit.y = init.y*fudge;
  2019.    }
  2020.  
  2021.    if(useinitorbit == 1)
  2022.       lold = linitorbit;
  2023.    else
  2024.       lold = linit;
  2025.  
  2026.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2027.    lold.y += lparm.y;
  2028.    return(1); /* 1st iteration has been done */
  2029. #else
  2030.    printf("Called long_mandel_per_pixel\n");
  2031.    exit(0);
  2032. #endif
  2033. }
  2034.  
  2035. int julia_per_pixel()
  2036. {
  2037.    /* julia */
  2038.  
  2039.    if(invert)
  2040.    {
  2041.       /* invert */
  2042.       invertz2(&old);
  2043.  
  2044.       /* watch out for overflow */
  2045.       if(bitshift <= 24)
  2046.      if (sqr(old.x)+sqr(old.y) >= 127)
  2047.      {
  2048.         old.x = 8;    /* value to bail out in one iteration */
  2049.         old.y = 8;
  2050.      }
  2051.       if(bitshift >  24)
  2052.      if (sqr(old.x)+sqr(old.y) >= 4.0)
  2053.      {
  2054.         old.x = 2;    /* value to bail out in one iteration */
  2055.         old.y = 2;
  2056.      }
  2057.  
  2058.       /* convert to fudged longs */
  2059.       lold.x = old.x*fudge;
  2060.       lold.y = old.y*fudge;
  2061.    }
  2062.    else
  2063.    {
  2064.       lold.x = lx0[col]+lx1[row];
  2065.       lold.y = ly0[row]+ly1[col];
  2066.    }
  2067.  
  2068.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2069.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2070.    ltmp = lold;
  2071.    return(0);
  2072. }
  2073.  
  2074. marks_mandelpwr_per_pixel()
  2075. {
  2076. #ifndef XFRACT
  2077.    mandel_per_pixel();
  2078.    ltmp = lold;
  2079.    ltmp.x -= fudge;
  2080.    LCMPLXpwr(lold,ltmp,ltmp);
  2081.    return(1);
  2082. #endif
  2083. }
  2084.  
  2085. int mandel_per_pixel()
  2086. {
  2087.    /* mandel */
  2088.  
  2089.    if(invert)
  2090.    {
  2091.       invertz2(&init);
  2092.  
  2093.       /* watch out for overflow */
  2094.       if(bitshift <= 24)
  2095.      if (sqr(init.x)+sqr(init.y) >= 127)
  2096.      {
  2097.         init.x = 8;  /* value to bail out in one iteration */
  2098.         init.y = 8;
  2099.      }
  2100.       if(bitshift >  24)
  2101.      if (sqr(init.x)+sqr(init.y) >= 4)
  2102.      {
  2103.         init.x = 2;  /* value to bail out in one iteration */
  2104.         init.y = 2;
  2105.      }
  2106.  
  2107.       /* convert to fudged longs */
  2108.       linit.x = init.x*fudge;
  2109.       linit.y = init.y*fudge;
  2110.    }
  2111.    else
  2112.       linit.x = lx0[col]+lx1[row];
  2113.    switch (fractype)
  2114.      {
  2115.     case MANDELLAMBDA:        /* Critical Value 0.5 + 0.0i  */
  2116.         lold.x = FgHalf;
  2117.         lold.y = 0;
  2118.         break;
  2119.     default:
  2120.         lold = linit;
  2121.         break;
  2122.       }
  2123.  
  2124.    /* alter init value */
  2125.    if(useinitorbit == 1)
  2126.       lold = linitorbit;
  2127.    else if(useinitorbit == 2)
  2128.       lold = linit;
  2129.  
  2130.    if(inside == -60 || inside == -61)
  2131.    {
  2132.       /* kludge to match "Beauty of Fractals" picture since we start
  2133.      Mandelbrot iteration with init rather than 0 */
  2134.       lold.x = lparm.x; /* initial pertubation of parameters set */
  2135.       lold.y = lparm.y;
  2136.       color = -1;
  2137.    }
  2138.    else
  2139.    {
  2140.       lold.x += lparm.x; /* initial pertubation of parameters set */
  2141.       lold.y += lparm.y;
  2142.    }
  2143.    ltmp = linit; /* for spider */
  2144.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2145.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2146.    return(1); /* 1st iteration has been done */
  2147. }
  2148.  
  2149.  
  2150. int marksmandel_per_pixel()
  2151. {
  2152.    /* marksmandel */
  2153.  
  2154.    if(invert)
  2155.    {
  2156.       invertz2(&init);
  2157.  
  2158.       /* watch out for overflow */
  2159.       if(sqr(init.x)+sqr(init.y) >= 127)
  2160.       {
  2161.      init.x = 8;  /* value to bail out in one iteration */
  2162.      init.y = 8;
  2163.       }
  2164.  
  2165.       /* convert to fudged longs */
  2166.       linit.x = init.x*fudge;
  2167.       linit.y = init.y*fudge;
  2168.    }
  2169.    else
  2170.       linit.x = lx0[col]+lx1[row];
  2171.  
  2172.    if(useinitorbit == 1)
  2173.       lold = linitorbit;
  2174.    else
  2175.       lold = linit;
  2176.  
  2177.    lold.x += lparm.x;     /* initial pertubation of parameters set */
  2178.    lold.y += lparm.y;
  2179.  
  2180.    if(c_exp > 3)
  2181.       lcpower(&lold,c_exp-1,&lcoefficient,bitshift);
  2182.    else if(c_exp == 3) {
  2183.       lcoefficient.x = multiply(lold.x, lold.x, bitshift)
  2184.      - multiply(lold.y, lold.y, bitshift);
  2185.       lcoefficient.y = multiply(lold.x, lold.y, bitshiftless1);
  2186.    }
  2187.    else if(c_exp == 2)
  2188.       lcoefficient = lold;
  2189.    else if(c_exp < 2) {
  2190.       lcoefficient.x = 1L << bitshift;
  2191.       lcoefficient.y = 0L;
  2192.    }
  2193.  
  2194.    ltempsqrx = multiply(lold.x, lold.x, bitshift);
  2195.    ltempsqry = multiply(lold.y, lold.y, bitshift);
  2196.    return(1); /* 1st iteration has been done */
  2197. }
  2198.  
  2199. marks_mandelpwrfp_per_pixel()
  2200. {
  2201.    mandelfp_per_pixel();
  2202.    tmp = old;
  2203.    tmp.x -= 1;
  2204.    CMPLXpwr(old,tmp,tmp);
  2205.    return(1);
  2206. }
  2207.  
  2208. int mandelfp_per_pixel()
  2209. {
  2210.    /* floating point mandelbrot */
  2211.    /* mandelfp */
  2212.  
  2213.    if(invert)
  2214.       invertz2(&init);
  2215.    else
  2216.       init.x = dx0[col]+dx1[row];
  2217.     switch (fractype)
  2218.       {
  2219.     case MAGNET2M:
  2220.         FloatPreCalcMagnet2();
  2221.     case MAGNET1M:         /* Crit Val Zero both, but neither   */
  2222.         old.x = old.y = 0.0; /* is of the form f(Z,C) = Z*g(Z)+C  */
  2223.         break;
  2224.     case MANDELLAMBDAFP:        /* Critical Value 0.5 + 0.0i  */
  2225.         old.x = 0.5;
  2226.         old.y = 0.0;
  2227.         break;
  2228.     default:
  2229.         old = init;
  2230.         break;
  2231.       }
  2232.  
  2233.    /* alter init value */
  2234.    if(useinitorbit == 1)
  2235.       old = initorbit;
  2236.    else if(useinitorbit == 2)
  2237.       old = init;
  2238.  
  2239.    if(inside == -60 || inside == -61)
  2240.    {
  2241.       /* kludge to match "Beauty of Fractals" picture since we start
  2242.      Mandelbrot iteration with init rather than 0 */
  2243.       old.x = parm.x; /* initial pertubation of parameters set */
  2244.       old.y = parm.y;
  2245.       color = -1;
  2246.    }
  2247.    else
  2248.    {
  2249.      old.x += parm.x;
  2250.      old.y += parm.y;
  2251.    }
  2252.    tmp = init; /* for spider */
  2253.    tempsqrx = sqr(old.x);  /* precalculated value for regular Mandelbrot */
  2254.    tempsqry = sqr(old.y);
  2255.    return(1); /* 1st iteration has been done */
  2256. }
  2257.  
  2258. int juliafp_per_pixel()
  2259. {
  2260.    /* floating point julia */
  2261.    /* juliafp */
  2262.    if(invert)
  2263.       invertz2(&old);
  2264.    else
  2265.    {
  2266.      old.x = dx0[col]+dx1[row];
  2267.      old.y = dy0[row]+dy1[col];
  2268.    }
  2269.    tempsqrx = sqr(old.x);  /* precalculated value for regular Julia */
  2270.    tempsqry = sqr(old.y);
  2271.    tmp = old;
  2272.    return(0);
  2273. }
  2274.  
  2275. int MPCjulia_per_pixel()
  2276. {
  2277. #ifndef XFRACT
  2278.    /* floating point julia */
  2279.    /* juliafp */
  2280.    if(invert)
  2281.       invertz2(&old);
  2282.    else
  2283.    {
  2284.      old.x = dx0[col]+dx1[row];
  2285.      old.y = dy0[row]+dy1[col];
  2286.    }
  2287.    mpcold.x = *pd2MP(old.x);
  2288.    mpcold.y = *pd2MP(old.y);
  2289.    return(0);
  2290. #endif
  2291. }
  2292.  
  2293. otherrichard8fp_per_pixel()
  2294. {
  2295.     othermandelfp_per_pixel();
  2296.     CMPLXtrig1(*floatparm,tmp);
  2297.     CMPLXmult(tmp,parm2,tmp);
  2298.     return(1);
  2299. }
  2300.  
  2301. int othermandelfp_per_pixel()
  2302. {
  2303.    if(invert)
  2304.       invertz2(&init);
  2305.    else
  2306.       init.x = dx0[col]+dx1[row];
  2307.  
  2308.    if(useinitorbit == 1)
  2309.       old = initorbit;
  2310.    else
  2311.       old = init;
  2312.  
  2313.    old.x += parm.x;     /* initial pertubation of parameters set */
  2314.    old.y += parm.y;
  2315.  
  2316.    return(1); /* 1st iteration has been done */
  2317. }
  2318.  
  2319. int MPCHalley_per_pixel()
  2320. {
  2321. #ifndef XFRACT
  2322.    /* MPC halley */
  2323.    if(invert)
  2324.       invertz2(&init);
  2325.    else
  2326.      init.x = dx0[col]+dx1[row];
  2327.  
  2328.    mpcold.x = *pd2MP(init.x);
  2329.    mpcold.y = *pd2MP(init.y);
  2330.  
  2331.    return(0);
  2332. #endif
  2333. }
  2334.  
  2335. int Halley_per_pixel()
  2336. {
  2337.    if(invert)
  2338.       invertz2(&init);
  2339.    else
  2340.       init.x = dx0[col]+dx1[row];
  2341.  
  2342.    old = init;
  2343.  
  2344.    return(0); /* 1st iteration is not done */
  2345. }
  2346.  
  2347.  
  2348. int otherjuliafp_per_pixel()
  2349. {
  2350.    if(invert)
  2351.       invertz2(&old);
  2352.    else
  2353.    {
  2354.       old.x = dx0[col]+dx1[row];
  2355.       old.y = dy0[row]+dy1[col];
  2356.    }
  2357.    return(0);
  2358. }
  2359.  
  2360. #if 0
  2361. #define Q0 .113
  2362. #define Q1 .01
  2363. #else
  2364. #define Q0 0
  2365. #define Q1 0
  2366. #endif
  2367.  
  2368. int quaternionjulfp_per_pixel()
  2369. {
  2370.    old.x = dx0[col]+dx1[row];
  2371.    old.y = dy0[row]+dy1[col];
  2372.    tmp.x = 0;
  2373.    tmp.y = 0;
  2374.    qc = param[0];
  2375.    qci = param[1];
  2376.    qcj = param[2];
  2377.    qck = param[3];
  2378.    return(0);
  2379. }
  2380.  
  2381. int quaternionfp_per_pixel()
  2382. {
  2383.    old.x = 0;
  2384.    tmp.x = 0;
  2385.    old.y = 0;
  2386.    tmp.y = 0;
  2387.    qc  = dx0[col]+dx1[row];
  2388.    qci = dy0[row]+dy1[col];
  2389.    qcj = param[2];
  2390.    qck = param[3];
  2391.    return(0);
  2392. }
  2393.  
  2394. int trigmandelfp_per_pixel()
  2395. {
  2396.    if(invert)
  2397.       invertz2(&init);
  2398.    else
  2399.       init.x = dx0[col]+dx1[row];
  2400.  
  2401.    if(useinitorbit == 1)
  2402.       old = initorbit;
  2403.    else
  2404.       old = init;
  2405.  
  2406.    old.x += parm.x;     /* initial pertubation of parameters set */
  2407.    old.y += parm.y;
  2408.    CMPLXtrig0(old,old);
  2409.    return(1); /* 1st iteration has been done */
  2410. }
  2411.  
  2412. int trigjuliafp_per_pixel()
  2413. {
  2414.    /* for tetrated types */
  2415.    if(invert)
  2416.       invertz2(&old);
  2417.    else
  2418.    {
  2419.       old.x = dx0[col]+dx1[row];
  2420.       old.y = dy0[row]+dy1[col];
  2421.    }
  2422.    CMPLXtrig0(old,old);
  2423.    return(0);
  2424. }
  2425.  
  2426. int trigXtrigmandelfp_per_pixel()
  2427. {
  2428.    if(invert)
  2429.       invertz2(&init);
  2430.    else
  2431.       init.x = dx0[col]+dx1[row];
  2432.  
  2433.    if(useinitorbit == 1)
  2434.       old = initorbit;
  2435.    else
  2436.       old = init;
  2437.  
  2438.    old.x += parm.x;     /* initial pertubation of parameters set */
  2439.    old.y += parm.y;
  2440.    CMPLXtrig0(old,tmp);
  2441.    CMPLXtrig1(old,tmp2);
  2442.    CMPLXmult(tmp,tmp2,old);
  2443.    return(1); /* 1st iteration has been done */
  2444. }
  2445.  
  2446. int trigXtrigjuliafp_per_pixel()
  2447. {
  2448.    if(invert)
  2449.       invertz2(&old);
  2450.    else
  2451.    {
  2452.       old.x = dx0[col]+dx1[row];
  2453.       old.y = dy0[row]+dy1[col];
  2454.    }
  2455.    CMPLXtrig0(old,tmp);
  2456.    CMPLXtrig1(old,tmp2);
  2457.    CMPLXmult(tmp,tmp2,old);
  2458.    return(0);
  2459. }
  2460.  
  2461. int MarksCplxMandperp(void)
  2462. {
  2463.    if(invert)
  2464.       invertz2(&init);
  2465.    else
  2466.       init.x = dx0[col]+dx1[row];
  2467.    old.x = init.x + parm.x; /* initial pertubation of parameters set */
  2468.    old.y = init.y + parm.y;
  2469.    tempsqrx = sqr(old.x);  /* precalculated value */
  2470.    tempsqry = sqr(old.y);
  2471.    Coefficient = ComplexPower(init, pwr);
  2472.    return(1);
  2473. }
  2474.  
  2475. QuaternionFPFractal()
  2476. {
  2477.    double a0,a1,a2,a3,n0,n1,n2,n3;
  2478.    a0 = old.x;
  2479.    a1 = old.y;
  2480.    a2 = tmp.x;
  2481.    a3 = tmp.y;
  2482.  
  2483.    n0 = a0*a0-a1*a1-a2*a2-a3*a3 + qc;
  2484.    n1 = 2*a0*a1 + qci;
  2485.    n2 = 2*a0*a2 + qcj;
  2486.    n3 = 2*a0*a3 + qck;
  2487.    /* Check bailout */
  2488.    magnitude = a0*a0+a1*a1+a2*a2+a3*a3;
  2489.    if (magnitude>rqlim) {
  2490.        return 1;
  2491.    }
  2492.    old.x = n0;
  2493.    old.y = n1;
  2494.    tmp.x = n2;
  2495.    tmp.y = n3;
  2496.    return(0);
  2497. }
  2498.  
  2499. /* -------------------------------------------------------------------- */
  2500. /*        Setup (once per fractal image) routines         */
  2501. /* -------------------------------------------------------------------- */
  2502.  
  2503. MandelSetup()        /* Mandelbrot Routine */
  2504. {
  2505.    if (debugflag != 90 && ! invert && decomp[0] == 0 && rqlim <= 4.0
  2506.        && bitshift == 29 && potflag == 0
  2507.        && biomorph == -1 && inside > -59 && outside >= -1
  2508.        && useinitorbit != 1)
  2509.       calctype = calcmand; /* the normal case - use CALCMAND */
  2510.    else
  2511.    {
  2512.       /* special case: use the main processing loop */
  2513.       calctype = StandardFractal;
  2514.       longparm = &linit;
  2515.    }
  2516.    return(1);
  2517. }
  2518.  
  2519. JuliaSetup()        /* Julia Routine */
  2520. {
  2521.    if (debugflag != 90 && ! invert && decomp[0] == 0 && rqlim <= 4.0
  2522.        && bitshift == 29 && potflag == 0
  2523.        && biomorph == -1 && inside > -59 && outside >= -1
  2524.        && !finattract)
  2525.       calctype = calcmand; /* the normal case - use CALCMAND */
  2526.    else
  2527.    {
  2528.       /* special case: use the main processing loop */
  2529.       calctype = StandardFractal;
  2530.       longparm = &lparm;
  2531.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2532.    }
  2533.    return(1);
  2534. }
  2535.  
  2536. NewtonSetup()        /* Newton/NewtBasin Routines */
  2537. {
  2538.    int i;
  2539.    extern int basin;
  2540.    extern int fpu;
  2541.    if (debugflag != 1010)
  2542.    {
  2543.       if(fpu != 0)
  2544.       {
  2545.      if(fractype == MPNEWTON)
  2546.         fractype = NEWTON;
  2547.      else if(fractype == MPNEWTBASIN)
  2548.         fractype = NEWTBASIN;
  2549.       }
  2550.       else
  2551.       {
  2552.      if(fractype == NEWTON)
  2553.            fractype = MPNEWTON;
  2554.      else if(fractype == NEWTBASIN)
  2555.            fractype = MPNEWTBASIN;
  2556.       }
  2557.       curfractalspecific = &fractalspecific[fractype];
  2558.    }
  2559.  
  2560.    /* set up table of roots of 1 along unit circle */
  2561.    degree = (int)parm.x;
  2562.    if(degree < 2)
  2563.       degree = 3;   /* defaults to 3, but 2 is possible */
  2564.    root = 1;
  2565.  
  2566.    /* precalculated values */
  2567.    roverd    = (double)root / (double)degree;
  2568.    d1overd    = (double)(degree - 1) / (double)degree;
  2569.    maxcolor    = 0;
  2570.    threshold    = .3*PI/degree; /* less than half distance between roots */
  2571. #ifndef XFRACT
  2572.    if (fractype == MPNEWTON || fractype == MPNEWTBASIN) {
  2573.       mproverd       = *pd2MP(roverd);
  2574.       mpd1overd    = *pd2MP(d1overd);
  2575.       mpthreshold  = *pd2MP(threshold);
  2576.       mpone       = *pd2MP(1.0);
  2577.    }
  2578. #endif
  2579.  
  2580.    floatmin = FLT_MIN;
  2581.    floatmax = FLT_MAX;
  2582.  
  2583.    basin = 0;
  2584.    if(roots != staticroots) {
  2585.       free(roots);
  2586.       roots = staticroots;
  2587.    }
  2588.  
  2589.    if (fractype==NEWTBASIN)
  2590.    {
  2591.       if(parm.y)
  2592.      basin = 2; /*stripes */
  2593.       else
  2594.      basin = 1;
  2595.       if(degree > 16)
  2596.       {
  2597.      if((roots=(CMPLX *)malloc(degree*sizeof(CMPLX)))==NULL)
  2598.      {
  2599.         roots = staticroots;
  2600.         degree = 16;
  2601.      }
  2602.       }
  2603.       else
  2604.      roots = staticroots;
  2605.  
  2606.       /* list of roots to discover where we converged for newtbasin */
  2607.       for(i=0;i<degree;i++)
  2608.       {
  2609.      roots[i].x = cos(i*twopi/(double)degree);
  2610.      roots[i].y = sin(i*twopi/(double)degree);
  2611.       }
  2612.    }
  2613. #ifndef XFRACT
  2614.    else if (fractype==MPNEWTBASIN)
  2615.    {
  2616.      if(parm.y)
  2617.      basin = 2; /*stripes */
  2618.       else
  2619.      basin = 1;
  2620.  
  2621.       if(degree > 16)
  2622.       {
  2623.      if((MPCroots=(struct MPC *)malloc(degree*sizeof(struct MPC)))==NULL)
  2624.      {
  2625.         MPCroots = (struct MPC *)staticroots;
  2626.         degree = 16;
  2627.      }
  2628.       }
  2629.       else
  2630.      MPCroots = (struct MPC *)staticroots;
  2631.  
  2632.       /* list of roots to discover where we converged for newtbasin */
  2633.       for(i=0;i<degree;i++)
  2634.       {
  2635.      MPCroots[i].x = *pd2MP(cos(i*twopi/(double)degree));
  2636.      MPCroots[i].y = *pd2MP(sin(i*twopi/(double)degree));
  2637.       }
  2638.    }
  2639. #endif
  2640.  
  2641.    param[0] = (double)degree; /* JCO 7/1/92 */
  2642.    if (degree%4 == 0)
  2643.       symmetry = XYAXIS;
  2644.    else
  2645.       symmetry = XAXIS;
  2646.  
  2647.    calctype=StandardFractal;
  2648. #ifndef XFRACT
  2649.    if (fractype == MPNEWTON || fractype == MPNEWTBASIN)
  2650.       setMPfunctions();
  2651. #endif
  2652.    return(1);
  2653. }
  2654.  
  2655.  
  2656. StandaloneSetup()
  2657. {
  2658.    timer(0,curfractalspecific->calctype);
  2659.    return(0);        /* effectively disable solid-guessing */
  2660. }
  2661.  
  2662. UnitySetup()
  2663. {
  2664.    periodicitycheck = 0;
  2665.    FgOne = (1L << bitshift);
  2666.    FgTwo = FgOne + FgOne;
  2667.    return(1);
  2668. }
  2669.  
  2670. MandelfpSetup()
  2671. {
  2672.    c_exp = param[2];
  2673.    pwr.x = param[2] - 1.0;
  2674.    pwr.y = param[3];
  2675.    floatparm = &init;
  2676.    switch (fractype)
  2677.    {
  2678.    case MANDELFP:
  2679.     /*
  2680.        floating point code could probably be altered to handle many of
  2681.        the situations that otherwise are using StandardFractal().
  2682.        calcmandfp() can currently handle invert, any rqlim, potflag
  2683.        zmag, epsilon cross, and all the current outside options
  2684.                              Wes Loewer 11/03/91
  2685.     */
  2686.     if (debugflag != 90
  2687.         && !distest
  2688.         && decomp[0] == 0
  2689.         && biomorph == -1
  2690.         && (inside >= -1 || inside == -59 || inside == -100)
  2691.         /* uncomment this next line if more outside options are added */
  2692.         /* && outside >= -5 */
  2693.         && useinitorbit != 1)
  2694.     {
  2695.        calctype = calcmandfp; /* the normal case - use calcmandfp */
  2696.        calcmandfpasmstart();
  2697.     }
  2698.     else
  2699.     {
  2700.        /* special case: use the main processing loop */
  2701.        calctype = StandardFractal;
  2702.     }
  2703.     break;
  2704.    case FPMANDELZPOWER:
  2705.       if(c_exp & 1) /* odd exponents */
  2706.          symmetry = XYAXIS_NOPARM;
  2707.       if(param[3] != 0 || (double)c_exp != param[2])
  2708.          symmetry = NOSYM;
  2709.       if(param[3] == 0.0 && debugflag != 6000 && (double)c_exp == param[2])
  2710.           fractalspecific[fractype].orbitcalc = floatZpowerFractal;
  2711.       else
  2712.           fractalspecific[fractype].orbitcalc = floatCmplxZpowerFractal;
  2713.       break;
  2714.    case MAGNET1M:
  2715.    case MAGNET2M:
  2716.       attr[0].x = 1.0;        /* 1.0 + 0.0i always attracts */
  2717.       attr[0].y = 0.0;        /* - both MAGNET1 and MAGNET2 */
  2718.       attrperiod[0] = 1;
  2719.       attractors = 1;
  2720.       break;
  2721.    case SPIDERFP:
  2722.       if(periodicitycheck==1) /* if not user set */
  2723.      periodicitycheck=4;
  2724.       break;
  2725.    case MANDELEXP:
  2726.       symmetry = XAXIS_NOPARM;
  2727.       break;
  2728. /* Added to account for symmetry in manfn+exp and manfn+zsqrd */
  2729. /*     JCO 2/29/92 */
  2730.    case FPMANTRIGPLUSEXP:
  2731.    case FPMANTRIGPLUSZSQRD:
  2732.      if(parm.y == 0.0)
  2733.         symmetry = XAXIS;
  2734.      else
  2735.         symmetry = NOSYM;
  2736.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  2737.         symmetry = NOSYM;
  2738.       break;
  2739.    case QUATFP:
  2740.       attractors = 0;
  2741.       periodicitycheck = 0;
  2742.       break;
  2743.    default:
  2744.       break;
  2745.    }
  2746.    return(1);
  2747. }
  2748.  
  2749. JuliafpSetup()
  2750. {
  2751.    c_exp = param[2];
  2752.    floatparm = &parm;
  2753.    if(fractype==COMPLEXMARKSJUL)
  2754.    {
  2755.       pwr.x = param[2] - 1.0;
  2756.       pwr.y = param[3];
  2757.       Coefficient = ComplexPower(*floatparm, pwr);
  2758.    }
  2759.    switch (fractype)
  2760.    {
  2761.    case JULIAFP:
  2762.     /*
  2763.        floating point code could probably be altered to handle many of
  2764.        the situations that otherwise are using StandardFractal().
  2765.        calcmandfp() can currently handle invert, any rqlim, potflag
  2766.        zmag, epsilon cross, and all the current outside options
  2767.                              Wes Loewer 11/03/91
  2768.     */
  2769.     if (debugflag != 90
  2770.         && !distest
  2771.         && decomp[0] == 0
  2772.         && biomorph == -1
  2773.         && (inside >= -1 || inside == -59 || inside == -100)
  2774.         /* uncomment this next line if more outside options are added */
  2775.         /* && outside >= -5 */
  2776.         && useinitorbit != 1
  2777.         && !finattract)
  2778.     {
  2779.        calctype = calcmandfp; /* the normal case - use calcmandfp */
  2780.        calcmandfpasmstart();
  2781.     }
  2782.     else
  2783.     {
  2784.        /* special case: use the main processing loop */
  2785.        calctype = StandardFractal;
  2786.        get_julia_attractor (0.0, 0.0);   /* another attractor? */
  2787.     }
  2788.     break;
  2789.    case FPJULIAZPOWER:
  2790.       if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2] )
  2791.          symmetry = NOSYM;
  2792.       if(param[3] == 0.0 && debugflag != 6000 && (double)c_exp == param[2])
  2793.           fractalspecific[fractype].orbitcalc = floatZpowerFractal;
  2794.       else
  2795.           fractalspecific[fractype].orbitcalc = floatCmplxZpowerFractal;
  2796.       get_julia_attractor (param[0], param[1]);    /* another attractor? */
  2797.       break;
  2798.    case MAGNET2J:
  2799.       FloatPreCalcMagnet2();
  2800.    case MAGNET1J:
  2801.       attr[0].x = 1.0;        /* 1.0 + 0.0i always attracts */
  2802.       attr[0].y = 0.0;        /* - both MAGNET1 and MAGNET2 */
  2803.       attrperiod[0] = 1;
  2804.       attractors = 1;
  2805.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2806.       break;
  2807.    case LAMBDAFP:
  2808.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2809.       get_julia_attractor (0.5, 0.0);    /* another attractor? */
  2810.       break;
  2811.    case LAMBDAEXP:
  2812.       if(parm.y == 0.0)
  2813.      symmetry=XAXIS;
  2814.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2815.       break;
  2816. /* Added to account for symmetry in julfn+exp and julfn+zsqrd */
  2817. /*     JCO 2/29/92 */
  2818.    case FPJULTRIGPLUSEXP:
  2819.    case FPJULTRIGPLUSZSQRD:
  2820.      if(parm.y == 0.0)
  2821.         symmetry = XAXIS;
  2822.      else
  2823.         symmetry = NOSYM;
  2824.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  2825.         symmetry = NOSYM;
  2826.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2827.       break;
  2828.    case QUATJULFP:
  2829.       attractors = 0;    /* attractors broken since code checks r,i not j,k */
  2830.       periodicitycheck = 0;
  2831.       break;
  2832.    default:
  2833.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2834.       break;
  2835.    }
  2836.    return(1);
  2837. }
  2838.  
  2839. MandellongSetup()
  2840. {
  2841.    FgHalf = fudge >> 1;
  2842.    c_exp = param[2];
  2843.    if(fractype==MARKSMANDEL && c_exp < 1)
  2844.       c_exp = 1;
  2845.    if(fractype==LMANDELZPOWER && c_exp < 1)
  2846.       c_exp = 1;
  2847.    if((fractype==MARKSMANDEL   && !(c_exp & 1)) ||
  2848.        (fractype==LMANDELZPOWER && c_exp & 1))
  2849.       symmetry = XYAXIS_NOPARM;    /* odd exponents */
  2850.    if((fractype==MARKSMANDEL && (c_exp & 1)) || fractype==LMANDELEXP)
  2851.       symmetry = XAXIS_NOPARM;
  2852.    if(fractype==SPIDER && periodicitycheck==1)
  2853.       periodicitycheck=4;
  2854.    longparm = &linit;
  2855.    if(fractype==LMANDELZPOWER)
  2856.    {
  2857.       if(param[3] == 0.0 && debugflag != 6000  && (double)c_exp == param[2])
  2858.           fractalspecific[fractype].orbitcalc = longZpowerFractal;
  2859.       else
  2860.           fractalspecific[fractype].orbitcalc = longCmplxZpowerFractal;
  2861.       if(param[3] != 0 || (double)c_exp != param[2] )
  2862.          symmetry = NOSYM;
  2863.     }
  2864. /* Added to account for symmetry in manfn+exp and manfn+zsqrd */
  2865. /*     JCO 2/29/92 */
  2866.    if((fractype==LMANTRIGPLUSEXP)||(fractype==LMANTRIGPLUSZSQRD))
  2867.    {
  2868.      if(parm.y == 0.0)
  2869.         symmetry = XAXIS;
  2870.      else
  2871.         symmetry = NOSYM;
  2872.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  2873.         symmetry = NOSYM;
  2874.    }
  2875.    return(1);
  2876. }
  2877.  
  2878. JulialongSetup()
  2879. {
  2880.    c_exp = param[2];
  2881.    longparm = &lparm;
  2882.    switch (fractype)
  2883.    {
  2884.    case LJULIAZPOWER:
  2885.       if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2])
  2886.          symmetry = NOSYM;
  2887.       else if(c_exp < 1)
  2888.          c_exp = 1;
  2889.       if(param[3] == 0.0 && debugflag != 6000 && (double)c_exp == param[2])
  2890.           fractalspecific[fractype].orbitcalc = longZpowerFractal;
  2891.       else
  2892.           fractalspecific[fractype].orbitcalc = longCmplxZpowerFractal;
  2893.       break;
  2894.    case LAMBDA:
  2895.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2896.       get_julia_attractor (0.5, 0.0);    /* another attractor? */
  2897.       break;
  2898.    case LLAMBDAEXP:
  2899.       if(lparm.y == 0)
  2900.      symmetry = XAXIS;
  2901.       break;
  2902. /* Added to account for symmetry in julfn+exp and julfn+zsqrd */
  2903. /*     JCO 2/29/92 */
  2904.    case LJULTRIGPLUSEXP:
  2905.    case LJULTRIGPLUSZSQRD:
  2906.      if(parm.y == 0.0)
  2907.         symmetry = XAXIS;
  2908.      else
  2909.         symmetry = NOSYM;
  2910.      if ((trigndx[0] == LOG) || (trigndx[0] == 14)) /* LOG or FLIP */
  2911.         symmetry = NOSYM;
  2912.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2913.       break;
  2914.    default:
  2915.       get_julia_attractor (0.0, 0.0);    /* another attractor? */
  2916.       break;
  2917.    }
  2918.    return(1);
  2919. }
  2920.  
  2921. TrigPlusSqrlongSetup()
  2922. {
  2923.    curfractalspecific->per_pixel =  julia_per_pixel;
  2924.    curfractalspecific->orbitcalc =  TrigPlusSqrFractal;
  2925.    if(lparm.x == fudge && lparm.y == 0L && lparm2.y == 0L && debugflag != 90)
  2926.    {
  2927.       if(lparm2.x == fudge)       /* Scott variant */
  2928.      curfractalspecific->orbitcalc =  ScottTrigPlusSqrFractal;
  2929.       else if(lparm2.x == -fudge)  /* Skinner variant */
  2930.      curfractalspecific->orbitcalc =  SkinnerTrigSubSqrFractal;
  2931.    }
  2932.    return(JulialongSetup());
  2933. }
  2934.  
  2935. TrigPlusSqrfpSetup()
  2936. {
  2937.    curfractalspecific->per_pixel =  juliafp_per_pixel;
  2938.    curfractalspecific->orbitcalc =  TrigPlusSqrfpFractal;
  2939.    if(parm.x == 1.0 && parm.y == 0.0 && parm2.y == 0.0 && debugflag != 90)
  2940.    {
  2941.       if(parm2.x == 1.0)    /* Scott variant */
  2942.      curfractalspecific->orbitcalc =  ScottTrigPlusSqrfpFractal;
  2943.       else if(parm2.x == -1.0)    /* Skinner variant */
  2944.      curfractalspecific->orbitcalc =  SkinnerTrigSubSqrfpFractal;
  2945.    }
  2946.    return(JuliafpSetup());
  2947. }
  2948.  
  2949. TrigPlusTriglongSetup()
  2950. {
  2951.    FnPlusFnSym();
  2952.    if(trigndx[1] == SQR)
  2953.       return(TrigPlusSqrlongSetup());
  2954.    curfractalspecific->per_pixel =  long_julia_per_pixel;
  2955.    curfractalspecific->orbitcalc =  TrigPlusTrigFractal;
  2956.    if(lparm.x == fudge && lparm.y == 0L && lparm2.y == 0L && debugflag != 90)
  2957.    {
  2958.       if(lparm2.x == fudge)       /* Scott variant */
  2959.      curfractalspecific->orbitcalc =  ScottTrigPlusTrigFractal;
  2960.       else if(lparm2.x == -fudge)  /* Skinner variant */
  2961.      curfractalspecific->orbitcalc =  SkinnerTrigSubTrigFractal;
  2962.    }
  2963.    return(JulialongSetup());
  2964. }
  2965.  
  2966. TrigPlusTrigfpSetup()
  2967. {
  2968.    FnPlusFnSym();
  2969.    if(trigndx[1] == SQR)
  2970.       return(TrigPlusSqrfpSetup());
  2971.    curfractalspecific->per_pixel =  otherjuliafp_per_pixel;
  2972.    curfractalspecific->orbitcalc =  TrigPlusTrigfpFractal;
  2973.    if(parm.x == 1.0 && parm.y == 0.0 && parm2.y == 0.0 && debugflag != 90)
  2974.    {
  2975.       if(parm2.x == 1.0)    /* Scott variant */
  2976.      curfractalspecific->orbitcalc =  ScottTrigPlusTrigfpFractal;
  2977.       else if(parm2.x == -1.0)    /* Skinner variant */
  2978.      curfractalspecific->orbitcalc =  SkinnerTrigSubTrigfpFractal;
  2979.    }
  2980.    return(JuliafpSetup());
  2981. }
  2982.  
  2983. FnPlusFnSym() /* set symmetry matrix for fn+fn type */
  2984. {
  2985.    static char far fnplusfn[7][7] =
  2986.    {/* fn2 ->sin     cos    sinh    cosh   exp    log    sqr  */
  2987.    /* fn1 */
  2988.    /* sin */ {PI_SYM,XAXIS, XYAXIS, XAXIS, XAXIS, XAXIS, XAXIS},
  2989.    /* cos */ {XAXIS, PI_SYM,XAXIS,  XYAXIS,XAXIS, XAXIS, XAXIS},
  2990.    /* sinh*/ {XYAXIS,XAXIS, XYAXIS, XAXIS, XAXIS, XAXIS, XAXIS},
  2991.    /* cosh*/ {XAXIS, XYAXIS,XAXIS,  XYAXIS,XAXIS, XAXIS, XAXIS},
  2992.    /* exp */ {XAXIS, XYAXIS,XAXIS,  XAXIS, XYAXIS,XAXIS, XAXIS},
  2993.    /* log */ {XAXIS, XAXIS, XAXIS,  XAXIS, XAXIS, XAXIS, XAXIS},
  2994.    /* sqr */ {XAXIS, XAXIS, XAXIS,  XAXIS, XAXIS, XAXIS, XYAXIS}
  2995.    };
  2996.    if(parm.y == 0.0 && parm2.y == 0.0)
  2997.     { if(trigndx[0] < 7 && trigndx[1] < 7)  /* bounds of array JCO 5/6/92*/
  2998.         symmetry = fnplusfn[trigndx[0]][trigndx[1]];  /* JCO 5/6/92 */
  2999.     }                 /* defaults to XAXIS symmetry JCO 5/6/92 */
  3000.    else
  3001.       symmetry = NOSYM;
  3002.    return(0);
  3003. }
  3004.  
  3005. LambdaTrigOrTrigSetup()
  3006. {
  3007. /* default symmetry is ORIGIN  JCO 2/29/92 (changed from PI_SYM) */
  3008.    if ((trigndx[0] == EXP) || (trigndx[1] == EXP))
  3009.       symmetry = XAXIS;
  3010.    if ((trigndx[0] == LOG) || (trigndx[1] == LOG))
  3011.       symmetry = XAXIS;
  3012.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3013.    return(1);
  3014. }
  3015.  
  3016. JuliaTrigOrTrigSetup()
  3017. {
  3018. /* default symmetry is XAXIS */
  3019.    if(parm.y != 0.0)
  3020.      symmetry = NOSYM;
  3021.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3022.    return(1);
  3023. }
  3024.  
  3025. ManlamTrigOrTrigSetup()
  3026. { /* psuedo */
  3027. /* default symmetry is XAXIS */
  3028.    if (trigndx[0] == SQR)
  3029.       symmetry = NOSYM;
  3030.    if ((trigndx[0] == LOG) || (trigndx[1] == LOG))
  3031.       symmetry = NOSYM;
  3032.    return(1);
  3033. }
  3034.  
  3035. MandelTrigOrTrigSetup()
  3036. {
  3037. /* default symmetry is XAXIS_NOPARM */
  3038.    if ((trigndx[0] == 14) || (trigndx[1] == 14)) /* FLIP  JCO 5/28/92 */
  3039.       symmetry = NOSYM;
  3040.    return(1);
  3041. }
  3042.  
  3043.  
  3044. ZXTrigPlusZSetup()
  3045. {
  3046. /*   static char far ZXTrigPlusZSym1[] = */
  3047.    /* fn1 ->  sin   cos    sinh  cosh exp   log   sqr */
  3048. /*         {XAXIS,XYAXIS,XAXIS,XYAXIS,XAXIS,NOSYM,XYAXIS}; */
  3049. /*   static char far ZXTrigPlusZSym2[] = */
  3050.    /* fn1 ->  sin   cos    sinh  cosh exp   log   sqr */
  3051. /*         {NOSYM,ORIGIN,NOSYM,ORIGIN,NOSYM,NOSYM,ORIGIN}; */
  3052.  
  3053.    if(param[1] == 0.0 && param[3] == 0.0)
  3054. /*      symmetry = ZXTrigPlusZSym1[trigndx[0]]; */
  3055.    switch(trigndx[0])
  3056.    {
  3057.       case COS:   /* changed to two case statments and made any added */
  3058.       case COSH:  /* functions default to XAXIS symmetry. JCO 5/25/92 */
  3059.       case SQR:
  3060.       case 9:   /* 'real' cos */
  3061.          symmetry = XYAXIS;
  3062.          break;
  3063.       case 14:   /* FLIP  JCO 2/29/92 */
  3064.          symmetry = YAXIS;
  3065.          break;
  3066.       case LOG:
  3067.          symmetry = NOSYM;
  3068.          break;
  3069.       default:
  3070.          symmetry = XAXIS;
  3071.          break;
  3072.       }
  3073.    else
  3074. /*      symmetry = ZXTrigPlusZSym2[trigndx[0]]; */
  3075.    switch(trigndx[0])
  3076.    {
  3077.       case COS:
  3078.       case COSH:
  3079.       case SQR:
  3080.       case 9:   /* 'real' cos */
  3081.          symmetry = ORIGIN;
  3082.          break;
  3083.       case 14:   /* FLIP  JCO 2/29/92 */
  3084.          symmetry = NOSYM;
  3085.          break;
  3086.       default:
  3087.          symmetry = XAXIS;
  3088.          break;
  3089.       }
  3090.    if(curfractalspecific->isinteger)
  3091.    {
  3092.       curfractalspecific->orbitcalc =  ZXTrigPlusZFractal;
  3093.       if(lparm.x == fudge && lparm.y == 0L && lparm2.y == 0L && debugflag != 90)
  3094.       {
  3095.      if(lparm2.x == fudge)       /* Scott variant */
  3096.          curfractalspecific->orbitcalc =  ScottZXTrigPlusZFractal;
  3097.      else if(lparm2.x == -fudge)  /* Skinner variant */
  3098.          curfractalspecific->orbitcalc =  SkinnerZXTrigSubZFractal;
  3099.       }
  3100.       return(JulialongSetup());
  3101.    }
  3102.    else
  3103.    {
  3104.       curfractalspecific->orbitcalc =  ZXTrigPlusZfpFractal;
  3105.       if(parm.x == 1.0 && parm.y == 0.0 && parm2.y == 0.0 && debugflag != 90)
  3106.       {
  3107.      if(parm2.x == 1.0)    /* Scott variant */
  3108.          curfractalspecific->orbitcalc =  ScottZXTrigPlusZfpFractal;
  3109.      else if(parm2.x == -1.0)    /* Skinner variant */
  3110.          curfractalspecific->orbitcalc =  SkinnerZXTrigSubZfpFractal;
  3111.       }
  3112.    }
  3113.    return(JuliafpSetup());
  3114. }
  3115.  
  3116. LambdaTrigSetup()
  3117. {
  3118.    int isinteger;
  3119.    if((isinteger = curfractalspecific->isinteger))
  3120.       curfractalspecific->orbitcalc =  LambdaTrigFractal;
  3121.    else
  3122.       curfractalspecific->orbitcalc =  LambdaTrigfpFractal;
  3123.    switch(trigndx[0])
  3124.    {
  3125.    case SIN:
  3126.    case COS:
  3127.    case 9:   /* 'real' cos, added this and default for additional functions */
  3128.       symmetry = PI_SYM;
  3129.       if(isinteger)
  3130.      curfractalspecific->orbitcalc =  LambdaTrigFractal1;
  3131.       else
  3132.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal1;
  3133.       break;
  3134.    case SINH:
  3135.    case COSH:
  3136.       symmetry = ORIGIN;
  3137.       if(isinteger)
  3138.      curfractalspecific->orbitcalc =  LambdaTrigFractal2;
  3139.       else
  3140.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal2;
  3141.       break;
  3142.    case SQR:
  3143.       symmetry = ORIGIN;
  3144.       break;
  3145.    case EXP:
  3146.       if(isinteger)
  3147.      curfractalspecific->orbitcalc =  LongLambdaexponentFractal;
  3148.       else
  3149.      curfractalspecific->orbitcalc =  LambdaexponentFractal;
  3150.       symmetry = XAXIS;
  3151.       break;
  3152.    case LOG:
  3153.       symmetry = NOSYM;
  3154.       break;
  3155.    default:   /* default for additional functions */
  3156.       symmetry = ORIGIN;  /* JCO 5/8/92 */
  3157.       break;
  3158.    }
  3159.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3160.    if(isinteger)
  3161.       return(JulialongSetup());
  3162.    else
  3163.       return(JuliafpSetup());
  3164. }
  3165.  
  3166. JuliafnPlusZsqrdSetup()
  3167. {
  3168. /*   static char far fnpluszsqrd[] = */
  3169.    /* fn1 ->  sin   cos    sinh  cosh    sqr    exp   log  */
  3170.    /* sin    {NOSYM,ORIGIN,NOSYM,ORIGIN,ORIGIN,NOSYM,NOSYM}; */
  3171.  
  3172. /*   symmetry = fnpluszsqrd[trigndx[0]];   JCO  5/8/92 */
  3173.    switch(trigndx[0]) /* fix sqr symmetry & add additional functions */
  3174.    {
  3175.    case COS: /* cosxx */
  3176.    case COSH:
  3177.    case SQR:
  3178.    case 9:   /* 'real' cos */
  3179.    case 10:  /* tan */
  3180.    case 11:  /* tanh */
  3181.    symmetry = ORIGIN;
  3182.     /* default is for NOSYM symmetry */
  3183.    }
  3184.    if(curfractalspecific->isinteger)
  3185.       return(JulialongSetup());
  3186.    else
  3187.       return(JuliafpSetup());
  3188. }
  3189.  
  3190. SqrTrigSetup()
  3191. {
  3192. /*   static char far SqrTrigSym[] = */
  3193.    /* fn1 ->  sin    cos    sinh   cosh   sqr     exp   log  */
  3194. /*         {PI_SYM,PI_SYM,XYAXIS,XYAXIS,XYAXIS,XAXIS,XAXIS}; */
  3195. /*   symmetry = SqrTrigSym[trigndx[0]];      JCO  5/9/92 */
  3196.    switch(trigndx[0]) /* fix sqr symmetry & add additional functions */
  3197.    {
  3198.    case SIN:
  3199.    case COS: /* cosxx */
  3200.    case 9:   /* 'real' cos */
  3201.    symmetry = PI_SYM;
  3202.     /* default is for XAXIS symmetry */
  3203.    }
  3204.    if(curfractalspecific->isinteger)
  3205.       return(JulialongSetup());
  3206.    else
  3207.       return(JuliafpSetup());
  3208. }
  3209.  
  3210. FnXFnSetup()
  3211. {
  3212.    static char far fnxfn[7][7] =
  3213.    {/* fn2 ->sin     cos    sinh    cosh  exp    log    sqr */
  3214.    /* fn1 */
  3215.    /* sin */ {PI_SYM,YAXIS, XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3216.    /* cos */ {YAXIS, PI_SYM,XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3217.    /* sinh*/ {XYAXIS,XYAXIS,XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3218.    /* cosh*/ {XYAXIS,XYAXIS,XYAXIS,XYAXIS,XAXIS, NOSYM, XYAXIS},
  3219.    /* exp */ {XAXIS, XAXIS, XAXIS, XAXIS, XAXIS, NOSYM, XYAXIS},
  3220.    /* log */ {NOSYM, NOSYM, NOSYM, NOSYM, NOSYM, XAXIS, NOSYM},
  3221.    /* sqr */ {XYAXIS,XYAXIS,XYAXIS,XYAXIS,XYAXIS,NOSYM, XYAXIS},
  3222.    };
  3223.    /*
  3224.    if(trigndx[0]==EXP || trigndx[0]==LOG || trigndx[1]==EXP || trigndx[1]==LOG)
  3225.       symmetry = XAXIS;
  3226.    else if((trigndx[0]==SIN && trigndx[1]==SIN)||(trigndx[0]==COS && trigndx[1]==COS))
  3227.       symmetry = PI_SYM;
  3228.    else if((trigndx[0]==SIN && trigndx[1]==COS)||(trigndx[0]==COS && trigndx[1]==SIN))
  3229.       symmetry = YAXIS;
  3230.    else
  3231.       symmetry = XYAXIS;
  3232.    */
  3233.    if(trigndx[0] < 7 && trigndx[1] < 7)  /* bounds of array JCO 5/22/92*/
  3234.         symmetry = fnxfn[trigndx[0]][trigndx[1]];  /* JCO 5/22/92 */
  3235.                     /* defaults to XAXIS symmetry JCO 5/22/92 */
  3236.    else {  /* added to complete the symmetry JCO 5/22/92 */
  3237.       if (trigndx[0]==LOG || trigndx[1] ==LOG) symmetry = NOSYM;
  3238.       if (trigndx[0]==9 || trigndx[1] ==9) { /* 'real' cos */
  3239.          if (trigndx[0]==SIN || trigndx[1] ==SIN) symmetry = PI_SYM;
  3240.          if (trigndx[0]==COS || trigndx[1] ==COS) symmetry = PI_SYM;
  3241.       }
  3242.       if (trigndx[0]==9 && trigndx[1] ==9) symmetry = PI_SYM;
  3243.    }
  3244.    if(curfractalspecific->isinteger)
  3245.       return(JulialongSetup());
  3246.    else
  3247.       return(JuliafpSetup());
  3248. }
  3249.  
  3250. MandelTrigSetup()
  3251. {
  3252.    int isinteger;
  3253.    if((isinteger = curfractalspecific->isinteger))
  3254.       curfractalspecific->orbitcalc =  LambdaTrigFractal;
  3255.    else
  3256.       curfractalspecific->orbitcalc =  LambdaTrigfpFractal;
  3257.    symmetry = XYAXIS_NOPARM;
  3258.    switch(trigndx[0])
  3259.    {
  3260.    case SIN:
  3261.    case COS:
  3262.       if(isinteger)
  3263.      curfractalspecific->orbitcalc =  LambdaTrigFractal1;
  3264.       else
  3265.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal1;
  3266.       break;
  3267.    case SINH:
  3268.    case COSH:
  3269.       if(isinteger)
  3270.      curfractalspecific->orbitcalc =  LambdaTrigFractal2;
  3271.       else
  3272.      curfractalspecific->orbitcalc =  LambdaTrigfpFractal2;
  3273.       break;
  3274.    case EXP:
  3275.       symmetry = XAXIS_NOPARM;
  3276.       if(isinteger)
  3277.      curfractalspecific->orbitcalc =  LongLambdaexponentFractal;
  3278.       else
  3279.      curfractalspecific->orbitcalc =  LambdaexponentFractal;
  3280.       break;
  3281.    case LOG:
  3282.       symmetry = XAXIS_NOPARM;
  3283.       break;
  3284.    default:   /* added for additional functions, JCO 5/25/92 */
  3285.       symmetry = XYAXIS_NOPARM;
  3286.       break;
  3287.    }
  3288.    if(isinteger)
  3289.       return(MandellongSetup());
  3290.    else
  3291.       return(MandelfpSetup());
  3292. }
  3293.  
  3294. MarksJuliaSetup()
  3295. {
  3296.    c_exp = param[2];
  3297.    longparm = &lparm;
  3298.    lold = *longparm;
  3299.    if(c_exp > 2)
  3300.       lcpower(&lold,c_exp,&lcoefficient,bitshift);
  3301.    else if(c_exp == 2)
  3302.    {
  3303.       lcoefficient.x = multiply(lold.x,lold.x,bitshift) - multiply(lold.y,lold.y,bitshift);
  3304.       lcoefficient.y = multiply(lold.x,lold.y,bitshiftless1);
  3305.    }
  3306.    else if(c_exp < 2)
  3307.       lcoefficient = lold;
  3308.    get_julia_attractor (0.0, 0.0);    /* an attractor? */
  3309.    return(1);
  3310. }
  3311.  
  3312. SierpinskiSetup()
  3313. {
  3314.    /* sierpinski */
  3315.    periodicitycheck = 0;        /* disable periodicity checks */
  3316.    ltmp.x = 1;
  3317.    ltmp.x = ltmp.x << bitshift; /* ltmp.x = 1 */
  3318.    ltmp.y = ltmp.x >> 1;            /* ltmp.y = .5 */
  3319.    return(1);
  3320. }
  3321.  
  3322. SierpinskiFPSetup()
  3323. {
  3324.    /* sierpinski */
  3325.    periodicitycheck = 0;        /* disable periodicity checks */
  3326.    tmp.x = 1;
  3327.    tmp.y = 0.5;
  3328.    return(1);
  3329. }
  3330.  
  3331. extern char usr_floatflag;
  3332.  
  3333. HalleySetup()
  3334. {
  3335.    /* Halley */
  3336.    periodicitycheck=0;
  3337.  
  3338.    if(usr_floatflag)
  3339.      fractype = HALLEY; /* float on */
  3340.    else
  3341.      fractype = MPHALLEY;
  3342.  
  3343.    curfractalspecific = &fractalspecific[fractype];
  3344.  
  3345.    degree = (int)parm.x;
  3346.    if(degree < 2)
  3347.       degree = 2;
  3348.    param[0] = (double)degree;
  3349.  
  3350. /*  precalculated values */
  3351.    AplusOne = degree + 1; /* a+1 */
  3352.    Ap1deg = AplusOne * degree;
  3353.  
  3354. #ifndef XFRACT
  3355.    if(fractype == MPHALLEY) {
  3356.       setMPfunctions();
  3357.       mpAplusOne = *pd2MP((double)AplusOne);
  3358.       mpAp1deg = *pd2MP((double)Ap1deg);
  3359.       mpdegree = *pd2MP((double)degree);
  3360.       mptmpparmy = *pd2MP(parm.y);
  3361.       mptmpparm2x = *pd2MP(parm2.x);
  3362.       mpone       = *pd2MP(1.0);
  3363.    }
  3364. #endif
  3365.  
  3366.    if(degree % 2)
  3367.      symmetry = XAXIS;   /* odd */
  3368.    else
  3369.      symmetry = XYAXIS; /* even */
  3370.    return(1);
  3371. }
  3372.  
  3373. StandardSetup()
  3374. {
  3375.    if(fractype==UNITYFP)
  3376.       periodicitycheck=0;
  3377.    return(1);
  3378. }
  3379.  
  3380.